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Abstract 


This  research  is  aimed  at  improving  the  state  of  the  art  of  GPS  algorithms,  namely, 
the  development  of  a  closed-form  positioning  algorithm  for  a  stand-alone  user  and  the 
development  of  a  novel  differential  GPS  algorithm  for  a  network  of  users. 

The  stand-alone  user  GPS  algorithm  is  a  direct,  closed-form,  and  efficient  new 
position  detennination  algorithm  that  exploits  the  closed-form  solution  of  the  GPS 
trilateration  equations  and  works  in  the  presence  of  pseudorange  measurement  noise  for 
an  arbitrary  number  of  satellites  in  view.  A  two-step  GPS  position  determination 
algorithm  is  derived  which  entails  the  solution  of  a  linear  regression  and  updates  the 
solution  based  on  one  nonlinear  measurement  equation.  In  this  algorithm,  only  two  or 
three  iterations  are  required  as  opposed  to  five  iterations  that  are  normally  required  in  the 
standard  Iterative  Least  Squares  (ILS)  algorithm  currently  used.  The  mathematically 
derived  stochastic  model-based  solution  algorithm  for  the  GPS  pseudorange  equations  is 
also  assessed  and  compared  to  the  conventional  ILS  algorithm.  Good  estimation 
performance  is  achieved,  even  under  high  Geometric  Dilution  of  Precision  (GDOP) 
conditions. 

The  novel  differential  GPS  algorithm  for  a  network  of  users  that  has  been  developed 
in  this  research  uses  a  Kinematic  Differential  Global  Positioning  System  (KDGPS) 
approach.  A  network  of  mobile  receivers  is  considered,  one  of  which  will  be  designated 
the  ‘reference  station’  which  will  have  known  position  and  velocity  information  at  the 
beginning  of  the  time  interval  being  examined.  The  measurement  situation  on  hand  is 
properly  modeled,  and  a  centralized  estimation  algorithm  processing  several  epochs  of 
data  is  developed.  The  effect  of  uncertainty  in  the  reference  receiver’s  position  and  the 


IV 


level  of  the  receiver  noise  are  investigated.  Monte  Carlo  simulations  are  performed  to 
examine  the  ability  of  the  algorithm  to  correctly  estimate  the  non-reference  mobile  users’ 
position  and  velocity  despite  substantial  satellite  clock  errors  and  receiver  measurement 
noise. 
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EFFICIENT  GPS  POSITIONING  DETERMINATION  ALGORITHMS 


I  Introduction 

The  NAVSTAR  Global  Positioning  System  (GPS)  is  a  space-based  satellite  radio 
navigation  system  that  provides  three-dimensional  (3-D)  user  positioning  by  solving  a  set 
of  nonlinear  trilateration  equations  using  pseudorange  measurements.  The  current 
method  of  solving  the  nonlinear  equations  is  to  linearize  the  pseudorange  equations  and 
calculate  the  user  position  iteratively,  starting  with  a  user-provided  initial  position  guess. 
In  this  research,  it  is  recognized  up  front  that  pseudorange  measurements  are  noise 
corrupted.  Hence,  the  stochastic  nature  of  the  measurements  is  reflected  in  the  GPS 
pseudorange  equations  from  the  onset  to  develop  a  probabilistically  sound  GPS  solution. 
By  stochastically  modeling  the  measurement  situation  at  hand,  solving  for  position  or 
velocity  becomes  a  stochastic  estimation  problem. 

This  research  consists  of  two  parts,  both  of  which  stochastically  model  the 
pseudorange  measurement  with  random  white  noise  and  solve  for  position  or  velocity  as 
a  stochastic  estimation  problem.  The  first  part  is  a  direct,  closed-form  and  efficient  new 
position  detennination  algorithm  that  exploits  the  closed-form  solution  of  the  GPS 
trilateration  equations  and  works  in  the  presence  of  pseudorange  measurement  noise  for 
an  arbitrary  number  of  satellites  in  view.  In  some  applications,  a  two-step  GPS  position 
determination  algorithm,  which  entails  the  solution  of  a  linear  regression  and  updates  the 
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solution  based  on  one  nonlinear  measurement  equation,  is  needed.  In  this  algorithm,  only 
two  or  three  iterations  are  required,  as  opposed  to  five  iterations  that  are  required  in  the 
standard  Iterative  Least  Squares  (ILS)  algorithm  currently  used.  The  mathematically 
derived  stochastic  model-based  solution  algorithm  for  the  GPS  pseudorange  equations  is 
also  assessed  and  compared  to  the  conventional  ILS  algorithm.  Good  estimation 
performance  is  achieved,  even  under  high  Geometric  Dilution  of  Precision  (GDOP) 
conditions. 

The  second  part  of  this  research  investigates  a  Kinematic  Differential  Global 
Positioning  System  (KDGPS)  algorithm.  A  number  of  mobile  receivers  is  considered, 
one  of  which  will  be  designated  the  ‘reference  station’,  which  will  have  known  position 
and  velocity  information  at  the  beginning  of  the  time  interval  examined.  The 
measurement  situation  on  hand  is  properly  modeled,  and  a  centralized  estimation 
algorithm  processing  several  epochs  of  data  is  developed.  The  effect  of  uncertainty  in  the 
reference  receiver’s  position  and  the  level  of  receiver  noise  are  investigated.  Monte 
Carlo  simulations  are  perfonned  to  correctly  estimate  the  non-reference  mobile  users’ 
position  and  velocity  despite  substantial  satellite  clock  errors  and  receiver  measurement 
noise. 

1. 1  Global  Position  ing  System 

The  Global  Positioning  System  consists  of  a  constellation  of  24  satellite  vehicles 
(SV)  arranged  in  six  orbital  planes  inclined  at  55  degrees  at  an  altitude  of  20,200km.  The 
constellation  continually  broadcasts  signals  that  can  be  utilized  by  a  receiver  on  the  user’s 
platform  based  on  the  concept  of  one-way  time  of  arrival  (TOA)  ranging.  The  GPS 
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receiver  determines  the  range  from  each  visible  satellite  to  the  user’s  platform.  This 
measured  range  is  called  a  “pseudorange”  since  there  are  errors  present  on  the  GPS 
signal. 


There  are  four  unknown  parameters  involved  with  GPS  positioning,  which 
include  3  Cartesian  position  parameters  x,y,  and  z,  and  GPS  clock  error.  To  determine  a 
solution  it  is  therefore  required  that  at  least  four  GPS  SVs  be  within  view  of  the  receiver. 
The  pseudoranges  from  the  SVs  are  used  to  determine  the  user’s  position  with  respect  to 
Earth.  A  typical  GPS  scenario  is  shown  in  Figure  1  below.  SV  geometry  plays  a  critical 
part  in  determining  GPS  positioning.  Poor  SV  geometry  with  respect  to  the  receiver 
produces  high  Geometric  Dilution  of  Precision  (GDOP),  which  can  adversely  affect  GPS 
position  solutions  [1]. 


Four  Unknowns:  X,  Y,  Z,  Coordinate  of  Receiver 
and  GPS  Provided  Time 


Figure  1.  A  GPS  Scenario  [2] 
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GPS  provides  two  types  of  services.  The  Standard  Positioning  Service  (SPS)  is 
designated  for  the  civilian  users.  The  Precise  Positioning  Service  (PPS)  is  intended  for 
U.S.  military  and  selected  government  agencies.  Access  to  PPS  signal  is  controlled 
through  two  cryptographic  features  denoted  as  Antispoofing  (AS)  and  Selective 
Availability  (SA).  AS  is  a  technique  intended  to  defeat  deception  jamming,  whereas  SA 
is  a  method  to  intentionally  inject  additional  error  onto  the  GPS  signal  to  deny  full  system 
accuracy. 

SA  decreased  the  positioning  accuracy  of  stand-alone  receivers  to  within  100- 
meters  RMS.  Military  receivers  utilized  de-encryption  techniques  to  remove  SA  and 
provide  position  accuracy  of  10-meters  root-mean-square  (RMS)  [1].  SA  was  effectively 
turned  off  in  2000. 

1.1.1  Stand-Alone  GPS  Positioning 

Single  point,  or  stand  alone,  GPS  techniques  utilize  signals  broadcast  from  the 
GPS  satellites  as  depicted  in  Figure  1.  First,  a  nominal  state  (consisting  of  user  position 
coordinates  and  receiver  clock  error)  can  be  represented  as 

A, 

yu 
K 

cStu 

At  a  given  measurement  epoch,  the  GPS  receiver  generates  a  set  of  n  pseudorange 
equations  (where  n  is  the  number  of  satellites  visible  to  the  receiver).  The  pseudorange 
from  the  user  receiver  to  the  ith  satellite  is  the  sum  of  the  true  range  plus  the  receiver 
range  equivalent  clock  error  (i.e.,  after  SV  clock,  ionospheric,  tropospheric,  etc.  errors 
have  been  corrected  or  deemed  negligible) 
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P;  =  v(x/  -xuf  +(y<  -yu)2  +(zi  ~zu)2  +cStu 


where 


(xu?  yu,  zu) 

(Xi,  yi,  Zi) 
c8tu 


=  ECEF  position  coordinates  of  the  user  (m) 

=  ECEF  position  coordinates  of  the  ith  satellite  (m) 
=  Range  equivalent  receiver  clock  error  (m) 


These  equations  are  non-linear,  and  several  techniques  have  been  developed  to 
solve  for  the  user  position.  These  include  closed-form  solutions,  Kalman  filtering,  and 
ILS  techniques  based  on  linearization  [1,  3].  Since  the  ILS  algorithm  is  arguably  the 
simplest  approach,  and  most  commonly  used,  it  will  be  subsequently  explained. 

Since  the  position  of  the  user  is  not  known,  an  estimate  of  the  user 
position  (jc,  y,  z)  is  used  to  generate  a  set  of  estimated  pseudoranges  to  each  of  the  n 
satellites,  i.e., 

Pi  =  V O ~i  -  xu  )2  +  ( V;  -  yu  )2  +  (z,  -  Zu  f  +  c8tu  (1-3) 

The  relationship  between  the  true  and  the  estimated  position  with  errors  can  be  written  as 

x„=x„+Ax„  (1-4) 

The  approximated  pseudorange  equations  (equation  (1-3))  are  then  linearized  using  a  first 
order  Taylor  series  approximation  to  yield 


Ap  =  HAx 


where 


AA  1  aA  Uyl  «zl  -  1 

Ap  =  :  H  =  :  !  i  ! 

AP„\  axn  ayn  -1 


with  the  elements  Apt  defined  as 
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z,  -  z 


AP>  =  Pi  -  Pi 


r,  =  V(*/  - *u )2  +  (y,-  - yJ2  +  (zi  - zu)2 


Solving  equation  (1-5)  has  the  solution 

Ax  =  H  'Ap  (1-7) 

The  values  obtained  for  Ax  are  used  to  update  equation  (1-4)  for  the  user  position. 

There  are  three  possible  cases  to  be  considered.  If  there  are  fewer  than  four 

satellite  pseudoranges  available,  the  position  cannot  be  determined  since  Ax  cannot  be 

resolved.  If  there  are  exactly  4  distinct  pseudoranges,  there  will  be  a  unique  solution. 
However,  if  there  are  more  than  four  satellites  visible,  as  is  generally  the  situation,  an 
overdetermined  linear  system  is  obtained,  and  no  solution  will  be  available  that  will 
perfectly  solve  the  equation  in  Ax.  For  this  case,  the  least  squares  solution  concept  can 
be  utilized. 

Basic  least-squares  technique  yields  the  solution 

Ax  =  (HrH)~'HrAp  (1-8) 

Alternatively,  a  weighted  least-squares  solution 

Ax  =  (HrCp1H)-1HrCp1AP  (1-9) 


can  be  used  when  the  pseudorange  measurements  have  different  error  statistics  or  when 
the  pseudorange  measurement  errors  are  correlated.  The  matrix  Cp  is  the  measurement 
error  covariance  matrix  (diagonal  terms  are  measurement  error  variances  and  off- 
diagonal  terms  are  correlation  between  measurement  errors).  It  is  noted  that  this 
weighted  solution  is  identical  to  the  unweighted  case  if  Cp  =  I  (identity  matrix). 
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For  the  over-determined  case,  there  is  generally  no  solution  for  Ax  that  exactly 
solves  the  measurement  equation.  However,  measurement  residuals,  v,  can  be  applied  to 
the  measurements  which  would  result  in 


Ap=HAx+u  (1-10) 

or 

v=Ap-HAx  (1-11) 

Single  point  positioning  estimates  only  receiver  clock  errors,  and  requires  a 
correction  for  the  satellite  clock  error.  Satellite  clock  error  corrections  can  be 


accomplished  as  described  in  ICD-GPS-200C  [3]: 


P corrected  ~  P  +cAt 


(1-12) 


where 


P  corrected 

P 

Atsv 

a®,  an,  ao,to 
Atr 
F 
e 

yfa 

Ek 


=  pseudorange  corrected  for  SV  clock  error 
=  original  (raw)  pseudorange  measurement 
=  SV  clock  correction  =  a®  +  afi  (t  -  t0)  +  a G(t  -  to)  +  Atr 
=  SV  clock  correction  parameters  from  navigation  message 

=  relativity  correction  =  Fe  Vu  sin  (Ek) 

=  constant  =  -4.442807633xl0'10  sec/(meter)1/2 
=  eccentricity  from  navigation  message 
=  square  root  of  semi-major  axis  from  navigation  message 
=  Eccentric  anomaly  (from  SV  position  calculation) 


Stand-alone  GPS  positioning  techniques  are  fast  and  reliable.  However,  the  poor 


accuracy  results  (typically  30m-50m)  are  undesirable.  In  order  to  obtain  the  relative 
accuracy  required  for  uses  such  as  vehicle  formation  flying,  kinematic  and  DGPS 


techniques  must  be  employed. 
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SV4 


SV2 


SV1 
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/cy 


SV3 


V  V  <S- 


Rjeference  Station 


Differential  Signal 


Mobile  User 


Figure  2.  A  Typical  DGPS  Situation 

1.1.2  Differential  GPS  Positioning  Techniques 

The  principal  idea  behind  DGPS  is  that  if  a  reference  receiver  is  available  with  a 
known  position,  common  errors  between  it  and  relatively  close  mobile  receivers  viewing 
the  same  satellites  can  be  removed.  Generally,  there  are  two  basic  DGPS  techniques, 
code  and  carrier-phase  based  [1]. 

1. 1. 3  Code-Based  Algorithm 

Let  the  reference  receiver,  m,  have  a  known  position  represented  as  (xm,  ym,  zm) 
and  the  reported  ith  satellite  position  (via  ephemeris  data)  be  represented  as  (x„  y,,  z,). 
The  geometric  distance  from  the  reference  receiver  to  the  i,h  satellite  is 

R‘m  =  VO,  ~xm)2  +(y,-ym)2  +(zi~zm)2  (1-13) 

The  reference  receiver  is  then  able  to  generate  a  pseudorange  measurement  to  the  ith 
satellite  as 

Pm  =  K,  +  £m,  space  +  8  m, user  +  C^m  (1-14) 
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where  s,„;Space  ,  £m, control  and  s„,,user  are  the  space,  control,  and  user  segment  induced 
pseudorange  errors,  respectively,  and  Stm  represents  the  reference  receiver’s  clock  offset 
from  GPS  system  time.  The  errors  are  summarized  in  Table  1  [1].  The  reference 
receiver  simply  resolves  the  difference  between  the  generated  pseudorange  to  the  ith 
satellite,  p‘m  and  it’s  geometric  rang e,R'n ,  to  create  the  differential  correction 


AP\n  =  P m  ~  K  =  £ 


m  m  m ,  space  °  m,  control  m,user 


+  £„ 


+  cSt„ 


(1-15) 


This  correction  term  is  utilized  by  the  user,  or  mobile,  receiver,  where  it  is 
differenced  with  the  users’  generated  pseudorange  measurement  to  the  same  satellite 


Pu  ~  AP‘m  =K+£U, space 


-(£ 


space 

+  £ 


+  £„ 


+  cdt„ 


m ,  space  m ,  control  m ,  user 


~\~  s  cSt  ) 

m . user  m  / 


(1-16) 


If  the  user’s  receiver  is  located  relatively  nearby  the  reference  receiver,  the  user’s 
receiver  pseudorange  equation  error  components  will  be  nearly  identical  to  those  of  the 
reference  receiver.  Exceptions  include  errors  that  are  not  common  to  both  receivers,  i.e. 
multi-path  and  receiver  noise.  Therefore,  the  corrected  user  pseudorange  is  obtained 

P u ,  corrected  ^ u  ^ u  ^ combined  ( 1  “  1  ^) 


where  s'u  is  the  residual  user  segment  error  (multi-path,  etc)  and  5tcombined  is  the  combined 
clock  offset  (dtu  -  Stm ) .  A  typical  comparison  between  stand-alone  measurements  and 
DGPS  measurements  is  given  in  Table  1  [1]. 
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Table  1.  Typical  Satellite  errors  before/after  DGPS  correction  [1] 


Segment  Source 

Error  Source 

GPS  la 
Error (m) 

Typical  GPS  la 
errors  after  DGPS 
corrections 

Space 

Satellite  clock  stability 

3.0 

0 

Satellite  perturbation 

1.0 

0 

Selective  Availability 

32.3 

0 

Other  (thermal,  radiation,  etc) 

0.5 

0 

Control 

Ephemeris  prediction  error 

4.2 

0 

Other  (thruster  performance,  etc) 

0.9 

0 

User 

Ionospheric  delay 

5.0 

0 

Tropospheric  delay 

1.5 

0 

Receiver  noise  and  resolution 

1.5 

2.1 

Multipath 

2.5 

2.5 

Other  (inter-channel  bias,  etc) 

0.5 

0.5 

System  UERE 

Total  (rms) 

33.3 

3.3 

1.1.4  Carrier  Phase-Based  Algorithm 

Obviously,  more  precise  position  information  can  be  obtained  by  measuring  the 
amount  of  shift  in  the  frequency  (Doppler  shift)  of  the  received  signal.  Typically,  this 
shift  in  carrier  frequency  arises  from  the  relative  motion  of  the  GPS  satellites  to  the  user 
resulting  in  Doppler  shift  frequencies  of  A f  =  ±  5000 Hz  with  respect  to  the  LI  and  L2 
carriers.  Thus, 

4 f  =  fR~fT  (1-18) 

where, 

/R  =  Frequency  received  at  the  receiver  (Hz) 
fj  =  Known  transmitted  frequency  (Hz) 

Integration  of  the  Doppler  shift  offset  over  time  can  result  in  extremely  precise 
measurements  (centimeter  range  for  the  LI  and  L2  frequencies).  Thus  the  carrier  phase 
measurements,  c can  be  calculated  by  integrating  the  Doppler  measurements  Afmeas 
over  the  time  epoch: 
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(1-19) 


0(t)  =  \¥meas(t)dt  +  <l>(to) 

h 

The  integer  portion  of  the  initial  carrier-phase  at  the  start  of  the  integration, 
is  referred  to  as  the  “carrier  phase  integer  ambiguity.”  This  integer  ambiguity  exists 
because  the  receiver  merely  begins  counting  carrier  cycles  once  the  user  tracks  the 
satellite  signal.  Resolution  of  this  integer  ambiguity  is  paramount  in  determining  the 
most  precise  range  measurement  possible.  Several  techniques  have  been  utilized  to 
resolve  this  problem,  most  popular  of  which  are  the  least-squares  iteration  process  or 
LAMBDA  methods  [4,  5,  6]. 

1.1.5  General  Kinematic  GPS  Techniques 


In  some  of  the  less  advanced  receivers,  user  velocity  is  calculated  as  the  time 
derivative  of  the  estimated  position,  i.e., 


du  u(t2)-u(tl) 
dt  t2  ( 


(1-20) 


In  general,  this  approach  yields  poor  results  and  is  acceptable  only  if  the  user’s  velocity  is 
constant  over  the  selected  time  interval. 


Many  receivers  process  Doppler  measurements  which  effectively  estimate  the 
Doppler  frequency  of  the  received  satellite.  The  satellite  velocity  vector  is  computed 
using  the  ephemeris  data  and  an  orbital  model  that  resides  within  the  receiver  [1].  At  the 
receiver  antenna,  the  received  frequency,/^,  is  given  by  the  classical  Doppler  equation 
(neglecting  relativistic  effects)  as  follows 

/»=/rd-4^L)  (1-21) 

c 
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where, 


fr=  transmitted  satellite  signal  frequency  (known) 
vr  =  satellite-to-user  relative  velocity  vector 

a  =  is  the  unit  vector  pointing  along  the  line  of  sight  from  the  user  to  the  satellite 
c  =  speed  of  light 

The  dot  product  represents  the  radial  component  of  the  satellite-user  relative  velocity 
vector  along  the  instantaneous  line  of  sight  to  the  satellite  vector,  a,  .  The  quantity  vr  is 
given  as  the  velocity  difference 

vr  =  v-u  (1-22) 

where  v  is  the  (known)  velocity  of  the  satellite  and  u  is  the  velocity  of  the  user  to  be 
determined  (both  referenced  to  a  common  ECEF  frame).  Therefore,  the  Doppler  offset 
due  to  the  relative  motion  satisfies 

4/'  =  A-/r=-/r(V~U)'a  U-23) 

c 

There  are  several  techniques  [1]  to  obtain  user  velocity,  u,  from  the  measured 
Doppler  frequency,  A f . 

For  the  jth  satellite,  equation  (1-23)  yields 

=/r,{  !-h(v, -«)•«,]}  (1-24) 

and  the  corrected  satellite  frequency  is  given  by 

fr  -  f.  '■  'V/  (1-25) 

where  fa  is  the  nominal  transmitted  frequency  and  A fT  is  the  correction  determined  from 
the  navigation  message  update. 
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The  measured  value  of  the  received  frequency  is  in  error  due  to  the  frequency  bias 
offset.  This  offset  is  related  to  the  drift  rate  of  the  user  clock,  relative  to  GPS  time,  by 

fRi=fj(\  +  Stu)  (1-26) 

where  8tu  is  considered  positive  if  the  user  clock  is  running  fast.  Through  algebraic 

manipulation,  this  can  be  rewritten  as 

c(Jj-fTj)  .  .  cfjStu 

- l-v  a  .  +  v  a  .  +  v  a  =x  a  +v  a  ■  +  z  a - -  ( 1  -27t 

/  xj  xj  yj  yj  zj  zj  u  xj  y  u  yj  u  zj  r  \A  ^  '  ) 

Tj  JTj 

where  Vj,  aj  is  the  jth  satellite  velocity  and  acceleration  component  respectively  and 
vu  =  (xu ,  yu ,  zu )  is  the  user  velocity.  To  simplify  this  equation,  we  let 


C(fj  ~  fjj  ) 


d.  = 

J  f v 


+  v  a  ■  +  v  a  .  +  v  a  . 

XJ  XJ  yj  yj  zj  zj 


(1-28) 


Since  the  tenn  / \  /  fT  is  =  1,  equation  (1-28)  can  be  written  as 


(1-29) 


d  ■  =  x  a  ■  +  v  a  ■  +  z  a  .  —  cSt 

j  u  xj  y  u  yj  u  zj  u 

We  now  have  four  unknowns  that  can  be  solved  for  by  using  measurements  from  four 
satellites  and  using  the  set  of  linear  equations 


d  = 


with  the  general  form 


which  can  be  solved  as 


d\ 

«vl 

«zl 

ll 

d2 

H  = 

«,2 

ay2 

«z2 

1 

g  = 

yu 

d3 

«,3 

ay3 

«z3 

1 

d4_ 

_«v4 

a  y4 

«z4 

1 

i 

1 

o 

ft?- 

52 

1 _ 

d  =  Hg 


(1-30) 


(1-31) 


g  =  H  d 


(1-32) 
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The  previously  stated  technique  for  obtaining  user  velocity,  equation  (1-32),  uses 
measurements  that  may  be  corrupted  by  measurement  noise  and/or  multi-path  errors.  A 
Kalman  Filter  method  may  be  used  to  compute  a  smoothed  navigation  solution. 

The  Kalman  filter  technique  is  a  recursive  algorithm  that  provides  optimal 
estimates  of  user  position,  velocity,  and  clock  drift  (PVT)  based  on  noise  statistics  and 
current  measurements.  The  filter  contains  a  dynamical  model  of  the  GPS  receiver 
platform  and  outputs  a  set  of  user  receiver  position  and  velocity  state  estimates  as  well  as 
the  associated  error  variances.  Kalman  Filters  entail  an  approach  which  simultaneously 
estimates  eight  states:  3  position  states,  3  velocity  states,  1  receiver  clock  bias  and  1 
receiver  clock  drift.  In  general,  the  velocity  estimates  are  only  valid  for  low  dynamic 
situations. 

Generally,  the  dynamical  model  can  be  derived  from  a  Taylor  series  expansion  of 
the  receiver  position  u  at  time  t  shortly  after  time  to: 


.  .  .  .  du(t)  1  du2(t)  2 

u{t)-u(t0)  +  — —  |,=,o  (t-t0)  +  —  |,=,o  (t-t0)  +••• 


dt 


2!  dt1 


1  du\t )  3 

...  -t - ; — (t  — t0)  +  h.o.t 

3!  dt 3  0 


(1-33) 


In  summary,  the  filter  propagates  the  platform  position  from  one  time  point  to  the 
next.  Using  these  propagated  states,  the  receiver  calculates  the  anticipated  pseudorange 
and  delta  pseudorange  (the  change  in  pseudorange  per  epoch  for  each  satellite).  Next,  the 
pseudorange  and  delta  pseudoranges  are  measured  and  the  difference  between  the 
anticipated  and  the  measured  values  (residuals  or  errors  in  the  user  position  and  velocity 
estimates)  is  taken.  These  errors  are  usually  sent  back  through  the  algorithm  to  be 
utilized  in  future  state  estimates. 
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Utilizing  a  Kalman  filter  allows  the  use  of  fewer  then  4  satellites  and  adjusts  the 
state  estimates  to  weight  the  effects  of  measurement  noise.  That  is,  when  measurement 
noise  is  high,  the  filter  places  heavier  weights  on  the  state  estimates  while,  on  the  other 
hand,  the  filter  places  heavier  weights  on  the  measurements  when  the  noise  is  low. 

1.2  Summary 

This  chapter  described  the  conventional  GPS  techniques  used  to  detennine  user 
position  and  velocity.  The  use  of  DGPS  techniques  allows  the  removal  of  nearly  all 
common  errors  in  each  of  the  three  segments  (User,  Control,  and  Space)  as  can  be  seen  in 
Table  1.  However,  errors  that  are  uncorrelated  from  receiver  to  receiver  are  not  removed 
and,  in  particular,  the  receiver  measurement  noise  is  actually  increased. 


15 


II  A  Direct ,  Closed-Form  and  Efficient  New  Position 

Determination  Algorithm 


This  Chapter  presents  the  development  of  the  closed-form  GPS  positioning 
algorithm  for  a  stand-alone  user.  First,  an  overview  of  the  closed-form  solution  is  given 
as  the  basis  for  the  underlying  positioning  concept.  Second,  the  theoretical  nature  of 
using  GPS  pseudorange  measurement  equations  in  the  presence  of  measurement  noise  for 
position  detennination  is  explored.  Third,  the  closed-fonn  solution  algorithm  is 
developed,  which  is  followed  by  the  development  of  a  Kalman-like  update  algorithm. 
The  closed-form  solution  for  a  scenario  with  four  satellites  in  view  is  then  examined. 
The  algorithm  is  then  summarized  step-by-step  outlining  how  the  Matlab  simulation  is 
developed.  This  chapter  ends  with  the  experiment  setup  and  simulation  results. 

2.1  Overview 

The  current  method  of  solving  for  GPS  user’s  position  is  to  linearize  the 
pseudorange  equations  and  calculate  the  user  position  iteratively,  starting  with  a  user- 
provided  initial  position  guess  [8].  For  near-earth  navigation,  the  center  of  the  earth  is  a 
good  initial  estimate,  and  the  currently  used  iterative  least  squares  (ILS)  algorithm 
converges  to  the  GPS  solution.  An  area  of  potential  improvement  that  has  been 
investigated  in  recent  years  is  the  use  of  non-iterative  closed-fonn  solutions  to  the 
nonlinear  pseudorange  GPS  equations.  Closed-fonn  solutions  have  been  developed  by 
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Bancroft  [9],  Leva  [10],  Krause  [11],  Abel  and  Chafee  [12,  13],  Hoshen  [14],  and  by 
Nardi  and  Pachter  [15,  16]. 

The  goal  of  this  research  is  to  develop  an  efficient  new  position  determination 
algorithm  that  uses  the  closed-fonn  solution  of  the  trilateration  equations  and  works  in 
the  presence  of  pseudorange  measurement  noise  for  an  arbitrary  number  of  satellites. 
The  new  position  determination  algorithm  will  focus  on  the  statistics  of  the  position 
estimates  and  employ  a  Kalman-like  filter. 

2.1.1  Closed-Form  Solutions 

The  closed-form  solution  developed  herein  is  an  improvement  over  [15]  and  [16] 
through  the  employment  of  a  more  rigorous  mathematical  fonnulation.  In  this  research 
and  in  our  previous  work  reported  in  [15]  and  [16],  an  over-determined  system  is  treated, 
making  use  of  all-in-view  (n  >  5)  satellites  as  opposed  to  using  just  four  satellites. 
Moreover,  this  work  departs  from  a  detenninistic  formulation  of  the  problem  ([9],  [11], 
[13],  and  [14])  and  specifically  addresses  the  development  of  a  reliable  closed-fonn 
solution  that  works  in  the  presence  of  measurement  noise.  Previous  works,  with  the 
exception  of  [17],  treated  the  pseudorange  equations  as  a  detenninistic  set  of  equations. 
In  [18],  the  deterministic  solution  of  [9]  is  adapted  to  account  for  measurement  noise,  and 
this  approach  is  further  developed  in  reference  [19]. 

In  this  dissertation,  pseudorange  measurements  are  recognized  to  be  noise 
conupted.  Hence,  the  stochastic  nature  of  the  measurements  is  reflected  in  the  GPS 
pseudorange  equations  from  the  onset  to  develop  a  probabilistically  sound  GPS  solution. 
By  stochastically  modeling  the  measurement  situation  at  hand,  solving  for  position 
becomes  a  stochastic  estimation  problem.  The  use  of  correct  stochastic  modeling  and 
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estimation  yields  a  GPS  solution  that,  in  addition  to  the  position  estimate,  provides  an 
estimate  of  the  measurement  noise  intensity,  provided  that  there  are  more  than  five 
satellites  in  view. 

Thus,  the  estimation  algorithm  developed  here  provides  a  data-driven  position 
(and  user  clock  bias)  estimation  error  covariance  prediction.  This  prediction  introduces  a 
new  confidence  factor  into  GPS  positioning  that  is  critical  for  the  downstream  integration 
of  GPS  and  Inertial  Navigation  Systems  (INS)  or  Synthetic  Aperture  Radar  (SAR) 
sensors.  Moreover,  an  attractive  feature  of  our  solution  is  its  good  estimation 
performance,  achieved  even  under  poor  GDOP  conditions  and  in  urban  environments 
where  the  number  of  visible  satellites  may  be  reduced  to  four. 

Moreover  a  direct,  or  autonomous,  solution  that  does  not  require  an  initial 
position  estimate  is  attractive  for  space  navigation  and  for  unusual  planar  array 
configurations  using  pseudolites,  where  the  iterative  process  is  sensitive  to  the  initial 
position  estimate  (e.g.,  the  application  discussed  in  references  [20]  and  [21]). 

Furthermore,  fast  solutions  that  require  fewer  iterations  and  Floating-Point  Operations 
(FLOPS)  are  attractive  for  high-speed  vehicles  such  as  spacecraft,  where  the 

computational  resources  may  be  limited. 

2.1.2  Development  of  the  Two-Step  Close-Form  Solutions 

The  method  of  linear  regression  from  statistics  is  used  to  obtain  preliminary 
closed-form  estimates  of  the  position  and  user  clock  bias.  The  number  of  in-view 

satellites  required  is  n  >  5.  In  addition,  a  data-driven  estimate  of  the  pseudorange 

measurement  noise  intensity  is  derived.  The  data-driven  estimation  of  the  measurement 
noise  intensity  requires  an  additional  satellite,  thus,  the  two-step  algorithm  developed  in 
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[22]  requires  at  least  six  satellites  in  view  (n  >  6).  In  section  2.4,  the  second  step  of  the 
new  algorithm  is  discussed  in  detail.  In  this  second  step,  the  closed-form  solution  is  used 
in  conjunction  with  one  nonlinear  measurement  equation;  thus,  an  update  step,  akin  to  a 
Kalman  filter  update  technique  is  developed.  This  supplementary  algorithm  uses  the 
solution  of  the  closed-form  algorithm  as  initialization.  The  two-step  algorithm  is 
validated  in  extensive  simulations.  Comparisons  are  drawn  with  results  achieved  using 
the  conventional  ILS  algorithm  currently  used  in  GPS  receivers.  Good  position  and  clock 
bias  estimates  are  obtained  using  the  two  step  algorithm  with  two  to  three  iterations  only, 
as  opposed  to  five  iterations  in  the  ILS  algorithm.  Also,  the  FLOPS  count  is  significantly 
reduced. 

2.2  Background 

This  section  is  concerned  with  the  theoretical  nature  of  a  position  determination 
algorithm  that  uses  the  closed-fonn  solution  of  the  trilateration  equations  and  works  in 
the  presence  of  pseudorange  measurement  noise  for  an  arbitrary  number  of  satellites. 
This  section  provides  the  basis  theory  that  underlines  the  problem. 

2.2.1  Pseudorange  Corruptions 

The  GPS  uses  the  radio  timing  principle  to  measure  range  between  the  satellites 
and  the  GPS  receiver,  making  it  a  time-of-arrival  system.  If  ranges  were  being  measured 
directly,  we  would  be  dealing  with  a  multilateration  system,  and  obtaining  a  position  fix 
would  be  easy.  Under  ideal  error-  and  noise-free  conditions,  if  both  the  satellite  and  the 
GPS  receiver’s  clock  were  perfectly  synchronized  with  GPS  time  with  no  error,  the 
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measured  range  would  be  the  true  range  [23].  However,  the  GPS  receiver  measures 
pseudoranges,  which  are  corrupted  by  the  receiver  clock  bias,  measurement  noise,  and 
other  error  sources.  The  latter  include  atmospheric  delays,  satellite  clock  errors, 
ephemeris  errors,  and  receiver- induced  errors. 

2. 2. 1.1  Receiver  Clock  Bias 

The  receiver  clock  bias  caused  by  the  difference  between  the  receiver  clock  time 
and  GPS  time  is  by  far  the  largest  contributor  to  the  difference  between  pseudorange  and 
range.  However,  the  receiver  clock  bias  is  common  to  the  set  of  simultaneous 
pseudorange  measurements,  enabling  it  to  be  treated  as  an  unknown  parameter  to  be 
estimated  along  with  the  user  position  coordinates;  hence  the  GPS  solution  consists  of  the 
user’s  three  space  coordinates  and  clock  bias. 

2.2. 1.2  Ephemeris  Corrections 

Ephemeris  corrections  provided  to  the  satellites  from  the  control  segment  could 
be  used  to  partially  eliminate  the  satellite  time  error  and  the  ephemeris  errors.  Known 
tropospheric  and  ionospheric  error  model  corrections  can  be  applied  to  partially 
compensate  for  tropospheric  and  ionospheric  delay  errors,  and  ionospheric  errors  can 
essentially  be  removed  using  dual-frequency  measurements.  Improved  receiver  designs 
are  used  to  minimize  the  effects  of  the  receiver-related  errors,  including  receiver  noise, 
code  loop  quantization  errors,  multipath  effects,  and  interchannel  errors. 

2.2.2  Pseudorange  Modeling  in  GPS 

If  the  residual  errors  are  grouped  together  under  one  random  variable,  v,  the  GPS 
pseudorange  equation  can  be  modeled  as  the  true  Euclidean  range  with  an  unknown  clock 
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bias  and  measurement  noise  superimposed;  thus,  the  stochastic  nonlinear  pseudorange 
measurement  equation  is 

R,  =  V (ux  -  *,•  )2  +  (u_v  ~  y,  f  +  k  -  z,  f  +b  +  v,  (2-1) 

This  equation  represents  the  ith  corrected  pseudorange  from  satellite  i,  i  =  1,  2,  3,  n, 
where  n  is  the  number  of  satellites  in  view;  (ux,  uy,  uz)  are  the  unknown  user  position  in 
earth-centered,  earth-fixed  (ECEF)  coordinates;  (x;,  y,,  z,)  are  the  known  coordinates  of 
the  ith  GPS  satellite  in  ECEF  coordinates;  b  is  the  unknown  range-equivalent  user  clock 
bias;  and  v,  is  zero-mean,  Gaussian,  pseudorange  measurement  noise.  It  is  reasonable  to 
assume  that  all  receiver  measurements  are  subject  to  the  same  receiver  noise  intensity; 
therefore,  they  will  have  the  same  variance,  a2.  However,  the  measurement  noise  terms 
Vi  are  not  correlated  between  satellites. 

Concerning  the  measurement  noise,  v:  Given  the  number  of  contributing  factors 
to  pseudorange  noise  and  the  lack  of  knowledge  of  their  characteristics,  it  is  reasonable  to 
assume  that  the  residual  pseudorange  noise  will  have  a  zero-mean  Gaussian  distribution 
by  invoking  the  Central  Limit  Theorem  which  states  that  the  sum  of  many  independent 
random  variables,  regardless  of  their  contribution,  will  approach  a  Gaussian  distribution 
[20].  The  Gaussian  pseudorange  noise  will  not  be  white  due  to  the  correlated  nature  of 
the  encompassed  errors  and  noise.  This  concern  is  alleviated  since  there  is  no 
requirement  for  the  pseudorange  measurements  to  be  uncorrelated  in  time,  because  the 
positioning  problem  will  be  treated  as  a  static  estimation  problem,  where  a  snap  shot  in 
time  is  treated  as  a  new  static  estimation  problem.  On  the  other  hand,  it  is  desirable  for 
the  development  of  the  stochastic  estimation  that  pseudorange  noise  be  uncorrelated 
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across  satellites.  This  assumption  will  not  hinder  the  development  of  the  stochastic 
closed-form  solution  to  the  GPS  pseudorange  equation  in  this  research,  since  the 
pseudorange  measurements  will  be  differenced,  thereby  eliminating  some  of  the  effects 
of  correlated  noise.  Thus,  the  pseudorange  measurement  noise  is  modeled  as 

v,-  e  jV(0,  cr2),  i  =  1,  2, ...,  n  (2-2) 

2.3  Summary  of  Results 

This  section  summarizes  the  theoretical  development  of  the  closed-form  solution 
for  user  position  estimation  using  GPS  pseudorange  measurements.  The  parameter 
vector  consisting  of  user  position  and  receiver  clock  bias  can  be  defined  as 

u  =(ux,uy,uz,b)T 

It  is  required  to  obtain  an  estimate  of  the  parameter  u  .  It  is  also  desirable  to  predict  the 
estimation  error  covariance.  The  GPS  position  calculation  algorithm  described  in 
Theorem  1  below  provides  the  steps  necessary  to  achieve  these  objectives.  The 
derivation  and  proof  of  correction  will  be  given  in  the  sections  that  follow. 

2.3.1  Theorem  1 

Assume  that  the  number  of  satellites  in  view  is  n  >  6.  Given  the  ephemeredes 
{(xi ,  v, ,  z.  )}"=1  of  the  n  satellites  and  the  n  pseudoranges  {ft  }"=1  data,  a  position  estimate 

and  a  prediction  of  its  estimation  error  covariance  can  be  obtained  in  two  steps: 

Step  1:  Form  the  (n- 1)  x  4  regressor  matrix 
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H  = 


0„  ~xx) 
(x„  -x2) 


( y  n  ~y\) 

(. yn-y2 ) 


(z„-z,) 

(z„-z2) 


i^-RJ 

(R2-K) 


•  •  •  • 

(xn  ~  Xn  I  )  (yn~y«- 1)  (Z«-Z»-l) 


the  measurement  vector  Z  ,  whose  elements  are 

Z,  =  | fe2  -  R,,2  +  x,,2  - xf  +  y2  - y2  +  z2  - z,.2 ),  i  =  1,  2, n  - 1, 
and  the  (n  -  1 )  x  (n  -  1)  covariance  matrix  R  1 ,  whose  elements  are 


O'-1) 


(*-')  ,r'i 


— ’  *  =  1,2,. 


Jfc=l 


n- 1 


>+Z' 


*•=1 


where 


r:  =  - 


Then  use  this  infonnation  to  obtain  the  preliminary  closed-form  parameter  estimate 


=  (htr1hYhtr1z 


Next  calculate  the  return  difference  vector 


Z  =  Z  -Hu 


The  pseudorange  measurement  noise  intensity  (a)  estimate  can  be  found  from 
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Ux0  U\  ’  UyO  U2  ’  U:d  l{3 
Form  the  (n- 1)  x  4  regressor  vector 


Then  form  a  scalar  measurement 

7  _  „  ,  (» VO  - *„ K  +  (« VO  -  )v„  +  Oh-0  -  K 

n  n  I  2  2  2* 

V(M-vO  -  Xn)  +  Ko  ~  y„  )  +  («z0  “  Z,X 
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and  the  (correlation)  vector  p 


p  =  cj2(r  - b  )(HTR  'H )  ‘H tR 


_l(n-l)jcl 


Then  the  intennediate  matrix  variable  Y  is 


Y  =  P7  +  ppT  +  — t—  (p:hpr  +  phTP: ) 

1-  p  n 


lhPTh) 


and  K  is  the  modified  Kalman  filter  gain,  given  by 


K  = 


1 


f 


1-  p  h 


1 


\ 


Yh-p 
\  +  h  Yh  j 


Update  according  to 


u+ =  u  +K(zn-hTu  ) 

p-;  =  {i-]k-pIhY+pVxv 


Iterate: 

Akin  to  the  iterated  Kalman  Filter  algorithm  used  in  Extended  Kalman  Filtering,  the 
measurement  Z„  and  the  vector  h  are  updated  about  the  improved  position  estimate.  The 
algorithm  iterates  using  the  preliminary  estimate  and  estimation  error  covariance 
available  prior  to  the  update  and  produced  by  the  linear  closed-form  algorithm.  Hence,  in 
the  h  and  Zn  formulae  only,  set 


Ux0  ~  Ux  ’  Uy0  -  Uy  ’  U zO  ~  11 
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and  repeat  the  calculation  of  u+and  P t+.  The  algorithm  is  hardwired  to  stop  after  3 

iterations  as  experiments  have  shown  that  the  solution  converges  in  2  or  3  iterations.  At 
the  conclusion  of  Step  2  the  information  on  the  parameter  extracted  from  the  data  is 

u  e  N  {u  + ,  P-+  ) 

The  derivation  and  proof  of  correctness  of  the  novel  GPS  positioning  algorithm  is 
outlined  in  the  next  sections,  where  the  respective  steps  1  and  2  are  discussed,  and  is  also 
adapted  to  the  case  where  only  four  satellites  are  in  view. 

2.4  Closed-Form  Solution 

The  derivation  of  the  algorithm  in  step  1  is  now  presented  in  detail.  Equation  (2- 
1)  can  be  written  as 

(«,  -  x,  )2  +  (uy  -yj  +  (u:  -  zt  )2  =  (R,  -b-Vf  f  (2-3) 

Expanding  equation  (3)  results  in  the  following 

2  2  2  2 

ur  +u„  +uv  -b  —2xiur—2yiuv—2ziu,+2Rib 

2  2  2  2  2  (2-4) 

=  R,  -  xi  -  y,  -  h  -  2Ry,  +  2bvi  +  v, 

It  is  noted  that  the  first  four  terms  in  equation  (2-4)  are  the  unknown  variables 
squared  and  that  they  are  common  to  all  n  equations.  This  presents  an  opportunity  for 
eliminating  the  nonlinear  terms  by  differencing;  hence,  the  nth  equation  is  subtracted  from 
the  remaining  (n  -  1)  equations.  The  resulting  (n  -  1)  equations  are  linear  in  the 
unknown  variables  and  can  be  expressed  as 
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(2-5) 


(*„  -  x,  )ux  +  (y„  -  yt  )uy  +  (z„  -  z,  )uz  +  (R,  -  R„  )b 

l/r>2  rt  2  ,  2  2.  2  2.  2  2) 

=  -(A  -  Rn  +  -  X,  +  yn  -  yt  +  z„  -  z,  j 

+  Rnvn  -  RtVi  -  bvn  +  bv,  + 1  (vf  -  v; ) 

As  a  by-product  of  the  preceding  operation,  the  nonlinear  nth  pseudorange 
equation  remains 

Rn  =  )2  +  (uy  -  y„  )2  +  («_-  -  )2  +  b  +  v„  (2-6) 

The  nr/'  equation  will  remain  unused  in  this  phase  of  the  development  but  will 
subsequently  be  used  as  an  additional  “measurement”  equation  in  the  Kalman  update  in 
step  2.  The  linear  regressions  in  equation  (2-5)  can  be  written  compactly  in  matrix 
notation  fonn  as 


Z  =  Hu +  V 


where  Z  is  the  measurement  vector,  given  by 


Z  = 


z 


n— 1 


and  its  element 


(2-7) 


Z 


xf+y„ 


2  ,  2 

■  V.  +  Z 

J  i  n 


H  is  the  (n—  1)  x  4  regressor  matrix  given  by 


(2-8) 
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H  = 


(*„-*  i)  (yn  -y,) 

(xn-x2)  (yn  -y2) 


(Z„~Z  l) 

(Zn-Zl) 


(Rt-RJ 

(Ri-RJ 


(xn~x„- 1)  (yn-yn- 1)  (z„-z„-i) 


and  u  is  the  vector  of  unknowns,  i7  =  [n  v ,  w  ,  wz ,  h]r  . 

Finally,  V  is  the  (n  -  1)  error  vector  given  by  V  =  [VI,  V2,  •  Vn_x ]r ,  where 


Vi  =  Rnvn  ~  R,vi  ~  b{v„  -  v, )+  y  (v,2  -  v„2 )  /  =  1,  2, ...,  n  - 1 


(2-9) 


To  obtain  an  estimate  of  u  ,  the  statistics  of  the  equation  error  V  must  be  derived  from 
the  known  statistics  of  the  pseudorange  measurement  noise  v, .  According  to  equation 
(2-2),  the  following  holds: 


£{v,}  =  0 

E {v2  }=  cr2  for  each  i  =  1,  2, ...,  n  —  1 

E\y.Vj)  =  Q(i*  j) 


Thus,  the  statistics  of  the  error  vector  V  are  calculated  as  being 


e{v,}  =  o 

a*  +  a Xr„  -h  f  +|«  -/>)) 


for  i  =  j 


for  i  *  j 


The  covariance  matrix 


Pp  =  E[VVT  J  can  be  expressed  in  tri-diagonal  form  as 


(2-10) 
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di  c  •  •  c 
c  d2  c  •  • 

•  c  d2  •  • 

•  c  c  •  c 

C  •  •  •  dn  , 

where 

and  the  diagonal  elements 

dj  =  a2  +  (Rn  -  b )2  +  (i?,  -Z))_  for  i  =  1,  2, n  —  1 
The  linear  regression  in  equation  (2-7)  is  used  to  obtain  an  estimate  of  the 

unknown  parameters  u  .  The  aim  is  to  obtain  u  that  minimizes  the  estimate  error  as 
weighted  by  the  inverse  covariance  of  the  equation  error  vector.  The  minimum  variance 
parameter  estimate  is 

il  =  (htP-1h)1HtP-1Z  (2-12) 

It  is  desirable  to  find  a  closed-form  solution  for  fj  1  to  reduce  the  computation 
load  of  our  GPS  positioning  algorithm.  If  P-  is  expressed  as  Py  =  co  2R  ,  where 
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1  •  •  1 


1  •  1 


1  1^.1 


(2-13) 


1  1  •  • 


Pfl=—Rl 

cu 


(2-14) 


The  elements  of  the  diagonal  of  R  are  a  function  of  b,  the  clock  bias  error,  and 
a  ,  the  standard  deviation  of  the  measurement  noise.  For  implementation  purposes,  it  is 

desirable  to  remove  this  dependency  before  finding  a  solution  for  R  1 . 

Since  a2  « (Rn  -b)2 ,  and  for  most  positioning  applications,  b  «(Rj  -bf ,  the 

diagonal  elements  of  R  can  be  simplified  as  shown  in  equation  (2-15).  To  further 
strengthen  the  validity  of  the  assumptions  made  to  form  equation  (2-15),  Rn  can  be 
selected  as  the  largest  of  all  available  pseudoranges.  Thus, 


d,  .  R2 

—  «  1  +  — ^ T 


(2-15) 


After  several  algebraic  steps  and  applying  the  Matrix  Inversion  Lemma,  it  is  found  that 


(*_,L  =4 


1  2 

1  r  i 

-  r. - 

^2^  i  n- 1 

<J  C  , 


,  i  =  1,2,. 


1  +  X7 
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\  1  rr 

R-'kj=—2-ri - >  i*j,j  =  \,2,...,n-\  (2-16) 

<j  c  .  \- 

1  +  Z'l- 

k=l 

where 

/  \  2 

,=f* 

UJ 

A  remarkable  property  of  the  estimate  in  equation  (2-12)  is  the  fact  that  it  is  not 
dependent  on  a ,  the  pseudorange  measurement  noise  variance.  Equation  (2-14)  yields 

the  equation  error  covariance  P-  as  simply  R  premultiplied  by  a  scalar  quantity.  In 

equation  (2-12),  the  scalar  premultiplier  of  P-  will  cancel  out;  therefore,  the  minimum 

variance  parameter  estimate  (position  and  clock  bias)  in  equation  (2-12)  can  be  rewritten 
in  an  equivalent  form  as 

a  =  {htRxh)1  HtRxZ  (2-17) 


Equation  (2-17)  is  used  for  coding  the  MATLAB  [24]  algorithm.  It  must  be  noted 

that  there  are  no  large  matrix  inversions  associated  with  this  solution  since  R  1  has  been 
determined  analytically  and  can  be  coded  directly  into  the  algorithm.  The  only  inversion 

that  needs  to  be  performed  is  that  of  the  small  (4  x  4)  matrix  [ht  RlH^j,  which  can  easily 
be  hardwired  into  the  GPS  receiver’s  algorithm. 

Furthermore,  it  follows  from  equation  (2-17)  that  the  covariance  of  the  estimate’s 


error  is  given  by 


(2-18) 
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Unlike  the  solution  estimate,  the  covariance  P-  is  dependent  on  cr ;  hence,  cr  must  be 


known  or  estimated  to  estimate  the  error  covariance. 


Substituting  equation  (2-7)  into  equation  (2-17)  yields 


=  (htr1hYhtr'(hu  +  v) 


—  u  + 


(htrah) 


H  R  V 


Now  we  substitute  equation  (19)  into  the  return  difference  equation: 


(2-19) 


Z  =Z-HS  =  Hu  +  V-  H\u  +  [htR1hY  HT  R  ~xv) 

=  -h{htr 1h) 'htR-1)v 

H  R  H)  HR,  then  the  return  difference  can 

be  expressed  as 


Z  =  MV 


(2-20) 


We  now  define  the  weighted  return  difference 


Z  =  R;'Z 


(2-21) 


where  Re  1  is  obtained  from  the  Cholesky  decomposition  of  R:R  =  RL,ReT  .  To  avoid  a 
large  matrix  inversion,  as  the  size  of  Re  is  (n  -  1)  x  (n  -  1),  Re  1  is  expressed  as  a 


function  of  R  : 


r;1=r/r1 


We  now  calculate  the  scalar  quantity 


z'z  =  zt(r;])t  r;'z  =  ztr'z  =  vtmtr~1mv 


It  can  be  shown  that 
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and 


Thus, 


e{ztz}=  ca2Tr(M ) 
Tr(M )  =  n-  5 


e{ztz|=  ccr2{n  -  5) 


(»-5) 

y 


(2-23) 


(2-24) 


Rearranging  equation  (2-24)  results  in  a  quadratic  equation  in  <j\  Solving  this 
quadratic  equation  and  substituting 

zTz  =  ztr'z 


yields  the  following  data-driven  estimate  of  the  measurement  noise  intensity  a2 : 


(7  «  -\K  - 


7 , 2(Z  Z) 
n  —  5 


(2-25) 


In  conclusion,  the  derived  linear  regression  of  equation  (2-7),  which  consists  of  (n 
-  1)  equations,  requires  that  (n  —  1)  be  at  least  4  to  provide  an  initial  estimate  of  the  four 
parameters  in  u  .  This  implies  that  a  minimum  satellite  availability  of  5  is  required  to 
produce  the  solution  given  in  equation  (2-17).  At  least  1  extra  satellite  is  required  for  the 
data-driven  prediction  of  the  estimation  error  covariance  P-  .  Thus,  a  minimum  satellite 


availability  of  6  is  required  to  produce  an  initial  estimate  of  the  four  parameters  in  u  and 
the  estimation  error  covariance.  If  the  pseudorange  measurement  noise  intensity  is 
known,  then  5  satellites  in  view  suffice. 


The  solution,  i.e.,  the  preliminary  parameter  estimate  u ,  is  based  on  n  -  I 
equations  only,  although  n  measurements  are  available  initially.  This  indicates  that  n 
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equations  should  be  used  to  obtain  the  parameter  estimate.  The  second  step  of  the 
algorithm  developed  in  this  dissertation  addresses  this  issue  as  it  uses  the  remaining  n,h 
pseudorange  equation  (see  equation  (2-6))  as  a  measurement  update  in  a  Kalman 
filtering-like  update  equation. 

2.5  Kalman  Update  Algorithm 

The  proof  of  the  correctness  of  the  algorithm  in  step  2  is  now  given.  The  concept 
behind  the  Kalman  update  solution  approach  is  similar  to  that  of  a  conventional  Kalman 

filter.  The  closed-form  solution  above  provides  a  preliminary  GPS  solution  estimate  ( u  ) 
and  the  associated  error  covariance  matrix  ( P:).  We  now  use  the  rih  equation  to  update 

the  previous  estimate  in  the  same  way  this  would  be  accomplished  during  the  update 
cycle  of  an  extended  Kalman  filter. 

The  approach  that  is  used  begins  with  the  linearization  of  equation  (2-6)  about  a 
nominal  position  estimate.  The  linearized  equation  is  then  manipulated  into  the  standard 
linear  measurement  form  as  described  in  [25]  and  used  to  update  the  estimate. 

Since  the  measurement  of  equation  (2-6)  is  nonlinear,  it  may  be  necessary  for  the 
process  to  continue  in  an  iterative  manner  until  convergence  within  a  predefined 
tolerance  is  achieved.  Simulation  results  show  that  the  algorithm  converges  after  two  or 
three  iterations. 

ft  is  noteworthy  that  the  Kalman  update  algorithm  presented  in  this  section  differs 
from  the  basic  Kalman  filter  [26,  27]  in  that  the  “new”  measurement  used  to  update  the 
preliminary  estimate  is  correlated  with  the  preliminary  estimate.  The  conventional 
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Kalman  filter  update  equation  does  not  allow  for  correlation  between  the  new 
measurement  and  the  previous  estimate;  hence,  a  Kalman-like  update  equation  that 
accommodates  this  correlation  and  that  addresses  the  specific  measurement  situation  on 
hand  needs  to  be  derived. 

The  first  step  in  the  mathematical  development  of  the  Kalman  update  algorithm  is 
to  linearize  equation  (2-6)  about  a  nominal  user  position  (uxo,  uyo,  uzo)  by  performing  a 
Taylor  series  expansion  and  neglecting  second-  and  higher-order  terms.  Through 
equation  manipulation  and  rearranging  and  redefining  of  terms,  the  following  equation  in 
the  form  of  a  linear  scalar  measurement  model  is  obtained: 


Z„  =  h'u  +  v_ 


(2-26) 


where  Zn  represents  the  scalar  measurement  and  is  defined  as 

(u  n—x  )x  +  [u  n  —  v  )v  +  [u  r,—z  )z 

^  ^  n  /  n  V  yO  S  n  Js  n  V  zO  n)  n 

V (W  VO  -  xn  )2  +  («  vO  -  y„  )  2+  (W_-0  -  Zn  f 

and  A  is  a  (4  x  1)  regressor  vector  defined  as 


5-! 

O 

1 

X 

a 

_ 1 

V 

0  -X,,)2  +  (uy0  - y,f  +  (M,0  - z„)2 
( uyo-yn ) 

V 

'(“xO  -  XnY  +  (Mv0  -  T„)2  +  Ofio  -  Znf 

V 

KUx  0  ~  Xnf  +  OfvO  -  y,,)2  +  Oho  -  Znf 

1 

(2-27) 


Recall  that  u  is  the  vector  of  unknowns  \lix  ,  uv ,  uz ,  bf ,  and  vn  is  the  pseudorange 
measurement  noise  associated  with  the  nth  measurement,  where  vn  ~  /V (() ,  cr2).  The 
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initial  nominal  user  position  [ux0,  uy0,  w_0)  is  provided  by  the  estimate  of  the  closed-form 
algorithm  in  step  1 : 

Uxo  =  Uy0  ~  fiy  ’  W-0  =  ^ : 

Equation  (2-26)  is  in  the  desired  linear  measurement  model  form  that  can  be  used  to 
update  the  solution  obtained  from  the  preliminary  closed-form  algorithm  in  a  Kalman- 
like  update  step. 

Now,  Zn  is  actually  part  of  the  measurements  used  to  obtain  the  closed-form 
solution  and  not  a  new  measurement,  as  would  be  the  case  in  a  conventional  Kalman 
filter  application.  The  noise  in  the  new  measurement  and  the  previously  derived  position 
estimation  error  are  therefore  correlated.  This  is  a  violation  of  the  basic  assumptions 
used  in  the  derivation  of  the  conventional  Kalman  filter  update  equations.  To  derive  the 
new  Kalman-like  update  equation,  it  is  necessary  to  know  the  relationship  between  the 
noise  in  the  new  measurement  (v„)  and  the  preliminary  estimate.  This  preliminary 

estimate  is  provided  by  the  closed-form  algorithm. 

The  linear  regression  used  for  the  closed-fonn  algorithm  was  defined  in  equation 

(2-7),  and  the  statistics  of  the  noise  vector  V  were  derived  in  equation  (2-9).  The  closed- 
form  algorithm  produced  an  estimate  of  the  GPS  unknown  parameters,  u ,  defined  in 
equation  (2-17),  and  an  estimate  of  its  covariance  matrix  ( Ps  ),  defined  in  equation  (2-18). 

Using  this  knowledge  of  the  estimated  GPS  solution,  the  true  GPS  parameter  vector  can 
be  expressed  as  the  random  variable 

u  -  u  +  W  (2-28) 

where  W  ~  N{ 0,  P- ) .  The  correlation  of  interest  between  vn  and  W  can  be  defined  as 
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(2-29) 


pSEf/Vv„}=E{v,w} 

To  determine  the  relationship  between  W  and  vn,  the  linear  regression  in 

equation  (2-7)  is  multiplied  from  the  left  by  HT R  1 ,  yielding  the  expression 

HTRlZ  =  HTRlHu  +  HrR{V  (2-30) 

Equation  (2-30)  can  be  solved  for  u  to  obtain 

u  =  {htRxhYhtR1Z-{htR-1hYhtR1V  (2-31) 

The  first  tenn  on  the  right  hand  side  of  equation  (2-31)  is  recognized  from 
equation  (2-17)  as  u  ;  therefore,  an  expression  for  W  in  terms  of  V  is  obtained: 

W  =  {htR1hY  HTRlV  (2-32) 

Furthermore,  an  expression  that  represents  the  variance  between  any  single  elements  of 
V  and  vn  is  detennined  by  exploiting  the  noise  statistics  of  V  derived  previously: 

Efcv„}  =  (r\Rn  ~ b ),  /  =  1,  2, ...,  n- 1  (2-33) 

Equation  (2-33)  yields  the  following  covariance  matrix: 

V 

1 

E{vv,}=a2(R,-b)''  (2-34) 

-  (h-1)j:1 

Using  equations  (2-32)  and  (2-34),  an  expression  for  the  covariance  between  W 
and  vn  is  determined: 
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1 

1 

p  =  cj2(R„ -b){HTR  xHYHTR-x  [  (2-35) 

-  {n—\  )x\ 

Next,  an  augmented  linear  regression  is  formulated  by  combining  equations  (2-26)  and 
(2-28).  The  augmented  linear  regression  is  expressed  as 

Za=Hau  +  Va  (2-36) 

where  Za  is  the  (5x1)  augmented  “measurement”  vector  defined  as 


Ha  is  the  (5  x  4)  augmented  regressor  defined  as 


and  Va  is  the  (5x1)  augmented  “measurement  noise”  vector  defined  as 


In  the  derivation  that  follows,  in  order  to  distinguish  the  preliminary  estimate  u 
and  P-  as  produced  by  the  closed-form  algorithm  from  the  estimate  that  will  be  obtained 
through  the  Kalman  update  the  following  notation  is  used: 

•  it  and  P-  represent  the  estimate  and  the  estimation  error  covariance  prior  to 
the  update,  respectively. 
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•  i7+and  P?  represents  the  estimate  and  the  estimation  error  covariance 
following  the  update,  respectively. 

To  obtain  the  updated  estimates  from  the  augmented  linear  regression  in  equation 
(2-36)  it  is  necessary  to  derive  the  covariance  of  the  augmented  noise  vector  Va .  Since 


the  statistics  of  the  noise  components  in  Va  have  already  being  determined,  the  equation 
error  covariance  matrix,  Ra ,  is  given  by 


R.  = 


Pr 


(2-37) 


The  updated  GPS  minimum  variance  solution  estimate  and  the  associated  covariance  are 
then  given  by  the  expressions: 


r  =p;  hi  r;'z. 
p:  =(ht,r;'h.Y 


(2-38) 

(2-39) 


The  expressions  in  equations  (2-38)  and  (2-39)  are  sufficient  to  obtain  the 
required  updates,  but  it  is  desirable  to  manipulate  and  reduce  the  equations  into  the  more 
familiar  and  computationally  efficient  fonn  of  the  classical  Kalman  filter  update  equa¬ 
tions.  After  lengthy  manipulations  and  applying  the  Matrix  Inversion  Lemma,  the 
Kalman-like  update  equations  in  the  desired  form  are  obtained: 

u+  =  u~  +  k(zh  - hTu )  (2-40) 

Pr;  =  {l-{{l-pTh)K  +  p},T}Y  (2-41) 


where  the  intermediate  variable  Y  is  the  modified  pre-update  covariance  matrix  given  by 

li  P-  h  - 1 


Y  =  pr;  + 


h  P-h- 1  T  1  (  T  TD  \ 

-r - 7—777  pp  H - —  [P-  hp  +  ph  P  ) 

(l  -pTh)2  \-pThXu  11  ’ 


(2-42) 
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and  K  is  the  modified  Kalman  filter  gain  given  by 


K  = 


1 


1  -pTh 


Yh-  p 


\  +  h  Yh 


(2-43) 


The  parameter  estimate  update  in  equation  (2-38)  appears  to  be  identical  to  that  of  the 
classical  Kalman  filter  update  equations.  However,  this  is  not  the  case  since  the  Kalman 
filter  gain,  equation  (2-43),  is  not  the  same. 

Remarks:  Note  that  in  the  special  case  of  the  classical  Kalman  Filter  with  no  correlation, 
p  =  0  and  Y  =  P-  .  For  this  special  case  the  classical  Kalman  Filter  update  formulae  are 


indeed  recovered: 


and 


K  = 


1 


-P:h 


\  +  hT  P~h 
u+  =  u  +  k(zh  -  hTU~ ) 


P-:={l-KhT)p: 


Equations  (2-40)  to  (2-43)  are  used  in  the  MATLAB  implementation  of  the 
Kalman  update  algorithm.  The  Kalman  update  algorithm  is  intended  to  refine  the  GPS 
closed-form  solution  estimate  in  a  direct  and  non-recursive  manner.  However,  since  the 
measurement  in  equation  (2-7)  is  nonlinear,  it  is  necessary  for  the  process  to  continue  in 
an  iterative  manner  until  convergence  within  a  predefined  tolerance  is  achieved.  Now 
recall  that  the  new  measurement  used  by  the  Kalman  update  algorithm  is  actually  the  nth 
pseudorange  equation  in  equation  (2-7)  being  linearized  about  the  position  estimate 
produced  by  the  closed-form  algorithm  in  step  1.  This  implies  that  how  well  the 
linearization  fits  the  true  unknown  GPS  parameters  is  dependent  on  how  good  the 
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solution  produced  by  the  closed-fonn  algorithm  is  to  begin  with.  To  alleviate  this 
undesired  dependency,  after  the  Kalman  update  algorithm  has  been  applied  once  and 
produces  an  improved  solution  estimate,  equation  (2-7)  is  once  again  linearized  about  the 
improved  position  estimate,  producing  a  new  and  better  linear  measurement  equation. 
This  is  akin  to  the  iterated  Kalman  filter  algorithm  used  in  extended  Kalman  filtering. 

The  Kalman  update  algorithm  is  applied  a  second  time  using  the  preliminary 
estimate  and  estimation  error  covariance  available  prior  to  the  update  and  produced  by 
the  linear  closed-form  algorithm,  not  the  solution  obtained  as  a  result  of  the  previous 
application  of  the  Kalman  update.  Theoretically,  this  process  can  be  continued 
recursively  until  convergence  to  the  best  possible  solution  is  achieved;  however,  for  the 
scenario  that  was  used  in  the  experiment,  after  two  or  three  applications,  the  change  in 
the  solution  estimate  was  insignificant.  Hence,  the  algorithm  was  hardwired  to  perform 
three  iterations. 


2.6  Position  Determination  Algorithm  for  Four  Satellites  in 
View 


In  previous  sections,  it  was  required  that  n  >  5.  The  derivation  in  step  1  of  the 
algorithm  for  the  case  of  four  satellites  in  view  is  now  investigated.  Setting  n  =  4, 
equation  (2-5)  becomes: 

(*4  -  K  +  (v4  -  y,  h  V  +  (Z4  -  Z;  k  +  (Ri  ~R4  )b 


2  2,2  2,2  2  \ 
*4  “  xi  +  T4  _  T;  +z4  ~  zi  ) 


(2-44) 


R4V4  -  RiVi  -  hv4  +  bv,  +  J  (f  -  V4  ) 
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The  linear  regression  in  equation  (2-44)  can  be  compactly  written  in  matrix  notation  fonn 


as 


Z  =  Hu +  V 


where  Z  is  the  measurement  vector,  now  given  by 


Z  = 


Zx 

Z2 

Z  3 


and  the  element 


Z,  = 


-{r?-r 

O  V  1 


2,2  2,2  2,2 
4  +x4-xi  +  y 4  -yt  +zA-z 


f) 


H  is  the  (3  x  4)  regressor  matrix  given  by 


H 

1 

1 _ 

(v4  -  yx) 

(z4~zl) 

(R.-RJ 

H  = 

(x4-x2) 

(y4  -  y2) 

(z4-z2) 

(R2-R4) 

_(x4-x3) 

(v4  -  V3) 

(z4  -Z3) 

(R3-R4) 

(2-45) 


(2-46) 


(2-47) 


Finally,  V  is  the  ( 1  x  3)  error  vector  given  by  V  =  [Vt ,  V2 ,  V3  ]T ,  where 

Vt  =  R4V4  -  Rivi  +  b(vl  -  v4 )  +  ^  (v,2  -  v2 )  /  =  1,  2,  and  3  (2-48) 

We  now  have  three  linear  equations  in  four  unknowns,  and  a  stand-alone 
parameter  estimate  can  not  be  obtained.  The  nonlinear  fourth  pseudorange  equation  is 
yet  to  be  used: 

R4  =  j (ux  -X4)2  +  (uy -y4f  +(uz-z4 )2  +  b  +  v4  (2-49) 

As  in  step  2  in  the  Kalman  update  algorithm,  the  fourth  equation  can  be  linearized  about 
a  nominal  user  position  (uxo,  uyo,  uzo)  by  performing  a  Taylor  series  expansion  and 


42 


neglecting  second-  and  higher-order  terms.  Through  equation  manipulation,  the 
following  equation  in  the  form  of  a  linear  scalar  measurement  model  is  obtained: 


z4  -  h  u  +  v4 


(2-50) 


where  z4  represents  the  scalar  measurement  and  is  defined  as 

(« .vo  -  *4  )x4  +  tfvo  -  y4  )y4  +  {uz  0  -  z4  )z 


Z4  —  R4  +  ' 


V(M-vO  -  xj  +  (Mv0  -y4f  +  («r0  -  Z4f 


and  A  is  a  4  x  1  vector  defined  as 


h  = 


(w,0  -x4) 


-JOfio  x4)  (Wv0  ~y4)  +  (Wz0  2 4) 

_ Ko~  V4) _ 

-\/(M.v0  _  T4)  (My0  ~~  T4)  (Mz0  _  Z4) 
(“zQ-^) 

■)J(UxO  _  '^*4  )  (Wy0  _  T4)  (Wz0  _  Z4) 

1 


(2-51) 


(2-52) 


Recall  that  u  is  the  vector  of  unknowns,  [ux,  u  ,  u.,  bf ,  and  V4  is  the  pseudorange 


measurement  noise  associated  with  the  fourth  measurement,  where  v4  ~  n((),  cr2).  The 
initial  nominal  user  position  is  assumed  to  be  at  the  origin  of  the  ECEF  coordinates 

Uxo  =  0’  Uy0  =  0»  Mz0  =  0 

Combining  equations  (2-45)  and  (2-51),  we  now  have  an  augmented  measurement  model 
represented  by 


z 

~h~\ 

~v~ 

= 

u  + 

_Z4_ 

y] 

_V4_ 

(2-53) 


=  H  u  +  V 

aug  aug 
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The  statistics  of  the  equation  error  V  and  v4  are  the  same  as  previously  derived  in  the 
closed  form  for  an  arbitrary  number  of  satellites  in  view.  Thus,  the  expected  value 
between  these  two  errors  can  be  calculated  as 

E\  (>.,!  =  ■'f i 11 1 ’  -  Ryy,  -v,)  + 1 (vf  - v42 ) v4 J 

=  RtE{  v42}  -R,E{  v4v4)  - { v42 }  - £■{  v,v4 }) 

-£{v2})  U-54) 

=  R4a2  -  ba1  +  i  (e { v,J }  £ { v4 }  -  E {  v4! }) 

=  a2(R4  -  b ),  i  =  1,  2,  and  3 

Using  equations  (2-10)  and  (2-54),  the  (4  x  4)  covariance  matrix  P-  can  be  expressed  as 

dx  c  c  cx 

„  c  d^  c  c, 

P,  =  a  2  1 

Kug  c  c  d3  cx 

cx  c  j  cx  1 

where 

c  =  -  +  (R4-b )  2 
cx  =  (R4  -b) 

and  the  diagonal  elements 

dx  -  <j2  +  (R4  -  bf  +  (f?,  -  bf  for  i  =  1,  2,  and  3 
Thus,  the  linear  regression  in  equation  (2-53)  can  be  used  to  obtain  an  estimate  of  the 

unknown  parameters  u  .  The  aim  is  to  obtain  u  that  minimizes  the  estimation  error  as 
weighted  by  the  inverse  covariance  of  the  equation  error  vector. 

The  minimum  variance  parameter  estimate  is 


(2-55) 

4x4 
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d  =  \H . Tp:}  H . )  H . TPr'  z. 


(2-56) 


To  reduce  the  computation  load  of  our  GPS  positioning  algorithm,  P-  is  expressed  as 


P-  =  cct2R  ,  where 


d, 

—  1  1 


c  c 

1  1  it 


(2-57) 


1  1  1  1  / 
—  —  —  l/c 

_  c  c  c 


Then 


p = — r- R 


(2-58) 


Again,  since  a2  (( ( R4  -bf ,  and  for  most  positioning  applications, 

the  diagonal  elements  of  R  can  be  simplified  as  shown  below  in  equation  (2-59).  To 
further  strengthen  the  validity  of  the  assumptions  made  to  form  equation  (2-59),  R4  can 
be  selected  as  the  largest  of  all  available  pseudoranges.  Thus, 


c  Ri  ■ 

Cj  1 

—  «  —  =  c, 

C  R, 


(2-59) 


(2-60) 


and  l/c  can  be  approximated  as 


1  1  ~2 

-*T7sci 
c  R , 


(2-61) 
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Equation  (2-61)  shows  that  the  error  covariance  P-  is  simply  R  premultiplied  by 

a  scalar  quantity.  In  equation  (2-56)  the  scalar  premultiplier  of  Py  will  cancel  out; 

therefore,  the  minimum  variance  parameter  estimate  in  equation  (2-56)  can  be  rewritten 
as 


(2-62) 


Equation  (2-62)  is  used  for  coding  the  MATLAB  [24]  algorithm.  Given  that  equation  (2- 
49)  is  nonlinear,  a  few  iterations  will  be  required  in  order  to  convert  the  estimated 

position  u  to  within  a  predefined  tolerance. 


2. 7  Experimental  Setup 

The  closed-form  linear  regression  algorithm  developed  in  this  dissertation 
requires  at  least  six  pseudoranges  measurements  to  produce  a  stand-alone  GPS  solution 
and  a  prediction  of  the  position  estimation  error  covariance.  In  terms  of  satellite 
availability,  the  worst-case  scenario  occurs  at  latitudes  in  the  range  of  35°  -  55°,  where 
there  are  six  or  fewer  GPS  satellites  available  20  percent  of  the  time.  However,  more 
than  six  satellites  are  in  view  most  of  the  time.  Therefore,  for  the  case  of  n  =  4,  the  four 
pseudorange  measurements  are  randomly  selected  from  the  available  pseudorange 
measurements.  For  n  >  5,  the  two-step  algorithm  uses  all  n  available  pseudorange 
measurements  to  produce  the  GPS  solution.  Satellite  availability  is  not  dependent  on  user 
position  longitude;  hence,  a  fixed  user  position  in  the  35°  -  55°  latitude  range  over  the 
continental  United  States,  40°  N  latitude,  105°  W  longitude,  at  an  altitude  of  300  m,  was 
selected. 
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The  geographic  coordinates  are  converted  to  ECEF  coordinates  and  used  to 
generate  the  experimental  datasets  using  GPSoft’s  Satellite  Navigation  Toolbox  for 
MATLAB  [28].  The  Satellite  Navigation  Toolbox  is  used  to  generate  realistic  GPS 
satellite  position  data  from  which  true  ranges  can  be  calculated  between  all  in-view  GPS 
satellites  and  the  position  of  the  receiver.  After  adding  a  range  equivalent  user  clock  bias 
of  1000  m  to  all  the  ranges,  a  zero  mean  random  noise  of  preselected  standard  deviation  a 
=  100  m  is  superimposed  to  represent  the  Gaussian  measurement  noise. 

The  Satellite  Navigation  Toolbox  has  the  capability  of  simulating  realistic  noise- 
corrupted  pseudorange  measurements  which  can  be  applied  directly  to  the  GPS  position 
determination  algorithm,  as  would  be  the  case  in  a  real-world  scenario.  In  our 
experiments,  the  GPSoft  toolbox  was  used  for  simulating  just  the  GPS  satellite  ephemeris 
data,  and  the  simulated  pseudoranges  were  generated  as  described  above.  This  was  done 
for  the  following  reasons: 

•  It  provides  a  more  structured  dataset  for  analysis  of  the  algorithm,  since  only 
the  desired  effects  are  being  considered  and  the  amount  of  noise  corruption  on  the 
pseudorange  measurements  is  exactly  controlled;  and 

•  Since  the  pseudoranges  are  produced  starting  from  exactly  known  position 
coordinates,  comparisons  against  the  true  position  for  determining  the  algorithm’s 
accuracy  are  now  possible. 

In  addition  to  using  the  novel  two-step  GPS  positioning  algorithm  developed  in 
this  dissertation,  calculations  were  also  performed  using  the  conventional  ILS  algorithm 
to  provide  a  comparison  baseline.  Similar  to  our  algorithm,  in  the  simulation  an 
enhanced  ILS  algorithm  from  [16]  which  uses  all  n  available  pseudorange  measurements 
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to  obtain  the  GPS  solution  is  exercised.  Thus,  the  regressor,  or  H  matrix,  is  the 
conventional  matrix  of  direction  cosines  with  a  vector  of  l’s  populating  the  last  column. 
The  (n  x  4)  H  matrix  is  a  tall  matrix,  so  the  generalized  inverse  is  used,  resulting  in  a 
least-squares  solution. 

2.8  Simulation  Results 

The  novel  algorithm  developed  in  this  research  was  tested  against  the 
conventional  ILS  algorithm.  The  results  discussed  in  this  dissertation  are  the  cumulative 
representation  of  5,000  Monte  Carlo  runs.  It  was  found  experimentally  that  5,000  Monte 
Carlo  runs  are  enough  for  the  average  miss  distance  and  its  standard  deviation  to 
converge  for  the  algorithms  discussed  in  this  dissertation.  To  provide  a  fair  comparison 
of  the  results  from  each  approach,  the  Gaussian  pseudorange  noise  realization  for  each 
satellite  is  kept  the  same  between  both  algorithms  for  any  given  Monte  Carlo  run. 

The  estimation  results  as  a  function  of  satellite  availability  are  plotted  in  Figure  2 
and  tabulated  in  Table  2.  “Miss  dist.”  is  the  experimentally  determined  3-D  distance 
between  the  true  user  position  and  the  estimated  position  averaging  over  5,000  Monte 
Carlo  runs,  and  “std(miss)”  is  the  experimentally  detennined  standard  deviation  of  “miss 
dist.”  over  the  5,000  Monte  Carlo  runs.  The  predicted  standard  deviation  of  the  miss 
distance  is  gauged  according  to 

Predicted  std(miss)=  Jp -+  +  P-'  +  P' 

where  P- +  is  the  estimation  error  covariance  matrix  provided  by  the  estimation  algorithm. 
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All  “miss  dist.”  results  have  been  normalized  with  respect  to  the  measurement 
noise  standard  deviation  a  .  a  is  the  average  of  the  predicted  values  of  a ,  and  std(  a  ) 
is  the  standard  deviation  of  this  average.  Both  &  and  std(  6  )  have  also  been  normalized 
with  respect  to  a  .  The  number  of  iterations  and  FLOPS  are  the  experimentally  recorded 
number  of  iterations  and  FLOPS  required  to  produce  the  solution,  averaged  over  the 
5,000  Monte  Carlo  runs. 

2. 8. 1  Iterative  Least  Squares  Algorithm  Benchmark 

The  experimental  average  miss  distance  and  its  standard  deviation  produced  by 
the  ILS  algorithm  are  used  as  a  baseline  for  comparison  to  the  algorithm  presented  in  this 
dissertation.  The  average  non-dimensional  miss  distance  is  a  function  of  the  number  of 
satellites  in  view  and  it  ranged  from  1.44  to  2.95  for  this  algorithm.  The  experimentally 
obtained  non-dimensional  standard  deviation  of  the  miss  distance  is  relatively  small  and 
it  ranged  from  0.72  to  1.65.  The  average  non-dimensional  miss  distance  for  the  new 
algorithm  ranged  from  1.43  to  2.95,  which  is  comparable  to  the  results  obtained  by  using 
the  ILS  algorithm. 

Irrespective  of  satellite  availability,  it  took  the  ILS  algorithm  5  iterations  to 
converge  to  the  required  threshold  for  accuracy.  During  our  experiment,  FLOPS  counts 
ranged  between  3119  and  5503  for  the  ILS  algorithm  while  the  close-fonn  algorithm 
produced  FLOPS  counts  between  3061  and  5017.  The  variance  in  FLOPS  is  a  function 
of  the  size  of  H,  which  changes  as  a  function  of  satellite  availability;  of  course  the  “miss 
dist.”  decreases  as  the  availability  of  satellites  in  view  increases. 
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Table  2.  Average  Results  from  5,000  Monte  Carlos  Runs 


n  =  4 

n  =  6 

«  =  7 

n  =  8 

n  =  9 

ILS 

Alg. 

New 

Alg 

ILS 

Alg. 

New 

Alg 

ILS 

Alg. 

New 

Alg 

ILS 

Alg. 

New 

Alg 

ILS 

Alg. 

New 

Alg 

a  /  <7 

N/A 

0.81 

0.89 

0.93 

0.94 

std(  a)/  <7 

N/A 

0.60 

0.47 

0.40 

0.35 

Experiment 
RMS(miss)  / 

<T 

N/A 

2.24 

1.92 

1.61 

1.63 

Predicted 
std(miss)  /  <7 

N/A 

2.15 

1.97 

1.64 

1.72 

No.  of 
Iterations 

5 

5 

5 

3 

5 

3 

5 

3 

5 

3 

FLOPS 

3139 

3061 

4115 

3080 

4535 

3675 

5013 

4194 

5503 

5017 

2.8.2  Novel  Algorithm  Results 

From  a  performance  point  of  view,  the  novel  two-step  algorithm  produced  results 
comparable  to  the  baseline  ILS  results.  As  shown  in  Table  2,  the  experiment  RMS  miss 
distances  resulted  from  the  two-step  algorithm  are  similar  to  those  yielded  by  the 
predicted  standard  deviation  of  the  miss  distance.  There  is  an  apparent  trend  that  as  more 
satellites  became  available,  the  error  got  smaller,  with  the  exception  of  an  abnormally  at  n 
=  9. 

Furthennore,  the  two-step  algorithm  takes  only  three  iterations  to  produce  a 
position  estimate  and  a  prediction  of  the  estimation  error  covariance,  while  the 
conventional  ILS  algorithm  takes  five  iterations  to  produce  only  the  position  estimate. 
As  a  result,  the  FLOPS  count  for  the  two-step  algorithm  is  consistently  lower  than  the 
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FLOPS  count  for  the  ILS  algorithm.  Concerning  the  n  dependence  of  the  FLOPS  count, 
we  see  that  the  miss  distance  decreases  as  n  increases.  Very  good  results  are  obtained  for 
n  >  7.  However,  even  though  the  FLOPS  count  of  both  algorithms  are  proportional  to 
satellite  availability  (n),  the  FLOPS  count  of  the  two-step  algorithm  increases  at  a  faster 
rate.  This  is  due  to  the  estimation  of  cr  in  the  novel  algorithm,  which  requires  operations 
on  (n  -  1  )x  (n-  1)  matrices. 

Concerning  the  number  of  iterations  in  the  novel  algorithm,  the  novel  algorithm 
could  be  hard  wired  to  only  two  iterations  instead  of  three  to  lower  the  FLOPS  count 
even  further.  Experimental  results  show  that  if  the  number  of  iterations  is  reduced  to 
two,  the  average  miss  distance  and  the  estimation  error  covariance  remain  practically 
unchanged  for  n  >  7  and  slightly  change  for  n  =  6. 

2. 8.2.1  Unconventional  Geometries 

It  is  interesting  to  exercise  the  algorithm  in  unconventional  high-GDOP  scenarios 
where  conventional  iterative  algorithms  tend  to  have  difficulties.  This  type  of  scenario 
can  be  expected  in  test  range  applications  ([20]  and  [21]).  In  this  dissertation,  the 
experimental  test  environment  consists  of  a  simulated  ground  planar  array  of  36 
pseudolites  evenly  spaced  in  a  circular  pattern  with  a  10,000  meters  radius  and  one 
pseudolite  at  the  center.  This  pattern  was  selected  because  it  represents  the  best 
achievable  ground  array  for  the  conventional  ILS  algorithm  if  the  user  is  directly  above 
the  center  pseudolite  [21].  This  number  of  pseudolites  was  selected  to  achieve  satellite 
availability  levels  that  allow  for  evaluation  of  the  algorithm’s  estimation  capability. 
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In  the  simulations,  the  user  position  was  10,000  m  directly  above  the  center  of  the 
circular  pattern,  and  the  center  pseudolite  was  moved  away  from  the  center  to  vary  the 
geometry.  In  this  test  environment,  the  conventional  ILS  algorithm  produced  fairly  good 
estimates  of  the  four  GPS  parameters  with  a  pseudolite  directly  below  the  user;  however, 
the  estimates  quickly  degraded  as  the  center  pseudolite  was  moved  away  from  directly 
below  the  user,  and  it  failed  to  produce  a  solution  when  the  center  pseudolite  was  offset 
by  more  than  400  m. 


Figure  3.  Optimum  Ground-Based  Planar  Array 

Unlike  the  results  obtained  for  the  typical  near-earth  GPS  scenarios,  the 
estimation  of  the  2-D  ux,  and  uy  user  position  coordinates  was  extremely  good,  with  mean 
errors  similar  to  those  obtained  with  the  conventional  ILS  algorithm;  however,  the  user 
altitude  (uz)  mean  estimation  error  was  very  large,  ranging  from  2.9  x  105  meters  to  1.1  x 


52 


106m.  Given  the  extremely  low  estimation  errors  in  the  ux  and  uy  user  position  coordinate 
estimates,  it  appears  that  the  geometry  produced  by  pseudolite  ground  planar  arrays  is 
bias.  More  research  will  be  needed  to  fully  characterize  any  potential  usefulness  of  the 
closed-form  algorithm  for  this  application. 

Table  3.  Ground  Test  Results  from  5,000  Monte  Carlos  Runs 


Scenario 

Centered  Tx  Offset  =  0m 

Centered  Tx  Offset 
=  400m 

Close-Form 

Estimated  a 

95.56 

95.22 

Error  x 

-0.53 

0.12 

Error  y 

0.64 

-0.14 

Error  z 

293416 

1113025 

Miss  distance 

293416 

1113025 

Kalman  -like 

Update 

Error  x 

-0.52 

0.17 

Error  y 

0.44 

-0.90 

Error  z 

-10152 

-9872 

Miss  distance 

-10152 

-9872 

Iterative  Least 
Squared 

Error  x 

-0.57 

0.27 

Error  y 

0.68 

0.19 

Error  z 

0.39 

11.46 

Miss  distance 

0.97 

11.47 

Taking  the  estimates  produced  by  step  1  and  applying  the  Kalman  update 
algorithm  in  step  2  improved  the  estimate  of  uz.  However,  the  error  in  uz  was  still  too 
large  to  render  its  estimate  useful.  Moreover,  ux  and  uy  were  slightly  affected,  sometimes 
for  the  worse  and  other  times  for  the  better. 

Based  on  the  results,  the  Kalman  update  (step  2)  does  not  prove  very  useful  in  this 
ground  planar  array  test  environment  as  it  does  not  provide  any  significant  improvement 
over  the  closed-form  algorithm  (step  1).  It  is  noted  that  the  user  altitude  (uz)  estimation 
error  also  increased  from  0.39  to  11.46  for  the  ILS  algorithm.  Hence,  both  the  ILS  and 
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the  two-step  algorithm  do  not  yield  good  estimates  when  our  planar  arrays  of  pseudolites 
is  used  and  the  center  pseudolite  is  offset  by  400  m  from  the  center  of  the  array. 
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Ill  Stochastic  Model-Based  DGPS  Estimation 

Algorithm 


This  chapter  presents  the  development  of  the  kinematic  differential  GPS 
algorithm  for  a  network  of  users,  the  experiment  setup,  and  simulation  results  from  a 
5,000  Monte  Carlos  runs. 

3.1  Introduction 

A  novel  kinematic  differential  GPS  algorithm  is  presented  to  process  GPS 
pseudorange  measurements  from  multiple  receivers  obtained  during  multiple 
measurement  epochs.  Specifically,  the  accurate  relative  (and  absolute)  positioning  of  a 
team  or  formation  of  mobile  vehicles  is  considered,  and  a  general  navigation  web  concept 
is  advanced.  The  measurement  situation  is  modeled  under  a  stochastic  framework,  and, 
rather  than  differencing  and  double  differencing  as  is  currently  practiced  in  conventional 
DGPS,  the  common  errors  are  explicitly  acknowledged  and  a  centralized  estimation 
algorithm  is  derived.  In  addition,  the  predicted  covariance  of  the  different  receivers 
position  estimation  errors  is  obtained.  Extensive  simulations  were  perfonned  to  validate 
the  novel  kinematic  DGPS  algorithm.  The  results  were  compared  to  conventional  DGPS 
scenarios,  where  the  reference  stations  are  stationary  or  there  was  only  one  available 
observation  epoch.  The  positioning  accuracy  improvements  achieved  were  gauged 
against  the  perfonnance  of  conventional  DGPS. 
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3.1.1  Background 

There  are  a  large  number  of  applications  requiring  the  ability  to  obtain  real-time 
kinematic  positioning  of  high  altitude  formation- flying  vehicles  which  include  Low  Earth 
Orbit  (LEO)  military  satellite  fonnations  for  earth  observation/surveillance  [29]. 
Therefore,  formations  of  satellites  might  replace  large  complex  satellites  and  allow  lower 
overall  cost  and  a  higher  degree  of  redundancy  [30].  Research  is  being  conducted  to 
replace  or  augment  the  existing  inertial  navigation  systems  on  board  satellites  with  GPS 
receivers  to  allow  meeting  the  stringent  positioning  requirements  for  the  formation. 

While  the  majority  of  current  research  is  focused  on  utilizing  carrier  phase  GPS 
techniques  to  acquire  position  and  velocity  information  for  each  epoch  [30,  31,  32,  33, 
34],  this  study  will  examine  the  general  principals  and  trends  in  utilizing  pseudorange 
measurements  only  over  a  series  of  measurement  epochs. 

3.1.2  Problem  Statement 

GPS  positioning  accuracy  is  limited  by  measurement  errors  that  can  be  classified 
as  either  common  mode  or  non-common  mode.  Common  mode  errors  have  nearly 
identical  effects  on  all  receivers  operating  in  a  limited  geographic  area,  as  is  the  case  in 
formation  flight.  The  common  mode  errors  are  the  satellite  ephemerides  errors,  the 
satellite  clock  errors,  and  the  effects  of  the  ionosphere  and  troposphere.  Non-common 
mode  errors  are  distinct  even  for  two  receivers  with  small  antennae  separation  and  consist 
of  receiver  noise  and  resolution,  multipath  effects,  and  inter-channel  bias,  etc. 

The  DGPS  navigation  concept  is  based  on  the  realization  that  close-by  (<  50  km) 
GPS  receivers  are  equally  affected  by  the  common  errors.  Conventional  DGPS  uses  a 
reference  station  at  a  known  ECEF  position  to  determine  corrections  that  other  local,  and 


56 


potentially  mobile,  GPS  receivers  within  close  range  (<  50km)  of  the  reference  station 
can  use  to  reduce  the  effects  of  GPS  common  mode  errors.  Standard  DGPS  navigation 
exploits  the  known  position  of  a  reference  station  and  the  existence  of  a  communications 
channel  to  the  moving  vehicle  (the  “rover")  to  broadcast  corrections  to  the  GPS  receiver 
on  the  rover,  thus  improve  the  latter's  positioning  accuracy.  In  its  most  basic  form  the 
DGPS  methodology  entails  the  application  of  reference  station  broadcast  differential 
corrections  to  the  user  (rover)  measured  pseudoranges.  Alternatively,  the  positioning 
error  measured  by  the  reference  receiver  is  directly  subtracted  from  the  rover’s  position 
measurement.  Thus  DGPS  yields  a  stand-alone  and  improved  user  (rover)  position 
estimate  [35]. 

3. 1. 3  Summary  of  Current  Knowledge 

Using  GPS  for  relative  and  absolute  positioning  of  formation  flying  vehicles  is 
relatively  new  [30,  31,  32,  33],  The  majority  of  the  current  research  involves  utilizing  the 
Carrier  Phase  DGPS  techniques  for  tracking  and  controlling  vehicles.  Simulated  results 
indicate  that  <  2cm  rms  position  errors  are  possible  [31]. 

3.1.4  Scope 

The  positioning  error  in  DGPS  is  caused  by  receiver  noise  and  resolution, 
multipath,  and  inter-channel  bias,  etc.  Multipath  error  is  addressed  by  carefully  choosing 
the  antenna’s  location,  using  choke  ring  antennae,  and  applying  advanced  signal 
processing  techniques.  Hence  there  is  a  strong  incentive  to  develop  methodologies  for 
mitigating  the  effects  of  measurement  noise  and  residual  errors  in  DGPS.  Obviously,  an 
approach  which  relies  on  averaging  out  the  random  measurement  noise  and  residual 
errors  is  called  for.  In  Kinematic  GPS  (KGPS),  the  mitigation  of  the  random  residual 
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errors  and  measurement  noise  induced  effects  can  be  addressed  by  several  approaches. 
One  approach  is  the  application  of  a  kinematic  model  to  user  motion  and  clock  error 
during  the  (short)  measurement  interval.  The  other  possible  approach  is  the  use  of  a 
measurement  record  obtained  over  multiple  observation  epochs.  The  third  potential 
approach  is  the  centralized  processing  of  the  GPS  pseudoranges  taking  into  account  the 
underlying  temporal  dependence  of  the  kinematic  variables.  Hence,  improved  user 
position,  velocity,  and  possibly,  acceleration  estimates  are  obtained.  This  improvement  in 
navigation  performance  is  obtained  irrespective  of  whether  differential  corrections,  as 
provided  by  DGPS,  are  applied  to  the  raw  pseudorange  measurements.  Thus,  in  this 
research  KGPS  is  used  to  mean  not  just  a  mobile  user,  but  also  the  application  of  a 
dynamic  model  to  rover  motion. 

Hence,  provided  that  a  communications  channel  from  the  reference  station  to  the 
user  is  established  and  a  measurements  record  obtained  over  multiple  observation  epochs 
is  available,  one  is  strongly  motivated  to  combine  KGPS  and  DGPS.  In  this  way, 
improved  user  position,  velocity  and  acceleration  estimates  are  obtained  by  virtue  of 
employing  KGPS  while,  at  the  same  time,  the  ill  effects  of  measurement  noise  and 
residual  errors  are  mitigated  and  the  maximum  benefit  is  extracted  out  of  DGPS.  In  this 
dissertation  the  concept  of  synergistic  ally  employing  DGPS  and  KGPS  navigation  in  a 
monolithic  computational  algorithm  is  undertaken.  The  developed  algorithm  is  referred  to 
as  Kinematic  Differential  GPS  (KDGPS).  Specifically,  the  DGPS  scenario  with  one 
reference  station  at  a  surveyed  position  is  referred  to  as  conventional  KDGPS.  The  DGPS 
scenario  with  multiple  reference  stations  at  known  locations  is  referred  to  as  augmented, 
or  network,  KDGPS.  When  the  navigation  of  a  team  of  mobile  vehicles  (rovers)  is 
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considered,  such  that  now  the  “reference  stations"  are  mobile  and  their  positions,  along 
with  the  user's  position,  are  not  known  and  need  to  be  determined,  we  then  refer  to 
relative,  or  generalized  KDGPS.  The  latter  is  the  object  of  this  research. 

In  conventional  DGPS,  a  “team"  consisting  of  two  members  only  is  considered, 
viz.,  the  user  and  the  reference  station.  Also,  a  one-way  communications  channel  between 
the  reference  station  and  the  user  is  established.  In  conventional  DGPS  the  reference 
station  is  stationary  and  it  is  located  at  a  known  (surveyed)  position.  When  navigating 
near  an  array  of  stationary  reference  stations  which  surround  the  rover,  the  ionospheric 
corrections  transmitted  to  the  rover  are  calculated  by  interpolating  the  corrections 
computed  by  the  GPS  reference  stations  -  one  then  refers  to  network  DGPS  [36].  The 
corrections  to  the  pseudorange  measurement  errors  are  obtained  using  the  conventional 
method  of  differencing  and  double  differencing.  All  available  a-priori  infonnation 
should  be  used.  The  ionospheric  error  can  be  reduced  using  modeling  and  interpolation. 
However,  since  the  many  GPS  users  use  single  frequency  receivers,  the  residual 
ionospheric  errors  remain  a  concern.  Furthermore,  predicted  satellite  ephemerides  and 
predicted  satellite  clock  error  corrections,  as  provided  by  NASA’s  Internet-based  Global 
Differential  GPS  (IGDG)  system  which  uses  a  world  wide  network  of  stationary 
reference  stations,  should  be  applied  [37],  [38].  At  the  same  time,  as  it  is  desirable  to 
move  from  meter  level  to  centimeter  level  positioning  in  the  relentless  quest  for  higher 
accuracy,  residual  errors  after  DGPS  correction  cannot  be  neglected. 

In  this  research  we  are  mainly  interested  in  a  scenario  where  the  GPS  receivers 
are  on  vehicles,  and  multiple  rovers  in  close  proximity  must  be  accurately  positioned 
relative  to  each  other,  as  in  satellite  or  aircraft  formation  flight  control.  One  then  refers  to 
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relative  GPS  [39,  40,  41].  However,  rather  than  using  the  conventional  differencing  and 
double  differencing  method  [39,  40,  41],  in  this  research  a  stochastic  modeling 
framework  is  used  where  the  common  errors  are  explicitly  acknowledged  and  treated  as 
parameters  to  be  estimated.  The  use  of  multiple  GPS  receivers  which  are  connected  by  a 
two-way  data  link  and  centralized  processing  of  the  measurements  using  statistical 
methods  of  parameter  estimation  has  the  potential  to  significantly  improve  the  quality  of 
the  navigation  solution.  The  correct  treatment  of  the  common  GPS  errors  in  a  centralized 
computational  algorithm  is  a  very  important  concept.  Thus,  in  this  research  ILS,  linear 
regression,  and  the  Singular  Value  Decomposition  (SVD)  method  are  used,  and  a  more 
general  DGPS  paradigm  is  established  where  the  concept  of  a  reference  station  is  done 
away  with  and  a  team  of  m  users  is  considered.  All  the  team  members  are  treated  as 
equals,  and  the  designation  of  a  reference  station  is  not  needed.  Moreover,  the 
requirement  that  one  member  of  the  team  is  stationary  and  at  a  known  position  is  done 
away  with.  In  our  estimation  algorithm,  the  inclusion  of  prior  infonnation  on  the 
reference  stations’  positions  is  allowed  for,  as  in  network  adjustment  in  surveying  [42], 
and  the  theory  of  conventional  DGPS  is  recovered.  In  addition,  the  predicted  covariance 
of  the  rovers’  positions  estimation  error  is  obtained. 

The  pseudoranges  to  the  n  satellites  in  view  measured  by  all  the  receivers  are 
communicated  to  a  central  processor  and  are  operated  on  by  a  centralized  (optimal) 
estimation  algorithm,  where  the  common  errors  are  properly  accounted  for,  thus 
obtaining  improved  estimation  performance.  The  central  processor  then  communicates 
the  correct  position  to  each  of  the  team  members.  Concerning  communication  needs:  The 
central  processor  operates  on  all  the  m  team  members'  (in  x  n)  pseudoranges,  and  returns 
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the  estimated  positions,  velocities,  and  accelerations,  of  each  of  the  m  vehicles.  The  use 
of  an  architecture  which  employs  a  centralized  processor  reduces  the  number  of  required 
two-way  communication  channels  from  m(m  -  l)/2  to  m .  Of  course,  the  centralized 
processor  could  reside  in  one  of  the  fonnation's  vehicles,  i.e.,  in  the  “leader",  in  which 
case  only  (m  -  1)  two-way  communication  channels  are  required.  Evidently, 

m  -1  <  1/2  m(m  -1)  for  m  >  2  (3-1) 

Thus,  a  novel  navigation  web  concept  is  advanced.  Also,  the  navigation  web 
development  is  motivated  by  self-calibration  measurement  methods  from  physics  and 
surveying.  For  example,  in  [21]  the  self-calibration  concept  was  exploited  for  GDOP 
reduction  in  a  pseudolites-based  measurements  scenario.  In  geodesy,  [42],  the  positioning 
accuracy  of  a  network  of  surveyed  landmarks  is  enhanced  using  additional  angle  and 
elevation  measurements.  However,  the  common  errors  in  the  pseudorange  measurements 
are  not  explicitly  acknowledged  and  are  addressed  using  the  conventional  differencing 
and  double  differencing  method.  In  this  dissertation,  the  application  of  stochastic 
modeling  to  DGPS  positioning  is  pursued,  and  the  required  mathematical  development  is 
presented.  Indeed,  improved  relative  and  absolute  position  estimates  for  all  the  team 
members  are  obtained. 

3.1.5  Main  Contribution 

For  a  special  case  of  a  conventional  DGPS  scenario  and  in  a  deterministic  world, 
the  results  obtained  by  standard  double  differencing  will  be  the  same  as  those  provided  in 
this  dissertation  -  if  the  effects  of  nonlinearity  caused  by  linearization  in  the  commonly 
used  IFS  algorithm  are  not  considered,  and  without  an  estimate  of  the  position  error. 
After  all,  if  it  all  comes  down  to  a  solution  of  a  linear  deterministic  system,  then  whether 
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one  differences  to  eliminate  the  common  variables  using  either  Gauss  elimination  or  uses 
matrix  inversion,  the  result  is  the  same  -  provided  that  numerical  problems  are  avoided. 
The  latter  is  guaranteed  if  the  condition  number  of  the  measurement  matrix  is  sufficiently 
low,  which  is  the  case  in  GPS  when  the  GDOP  is  low. 

In  practice,  the  pseudorange  measurements  are  corrupted  by  random  measurement 
noise  and  residual  errors.  Hence,  in  this  research  the  emphasis  is  on  the  stochastic 
modeling  of  the  DGPS  measurement  situation  at  hand.  By  falling  back  on  the  classical 
statistical  theory  of  linear  regression,  ILS  and  SVD  method,  the  DGPS  problem  is 
correctly  fonnulated  and  a  state  (e.g.  position)  estimate  is  rigorously  obtained,  including 
a  prediction  of  the  variance  of  the  estimation  error.  Moreover  the  SVD  method  adapted 
from  numerical  analysis  yields  insight  into  the  structure  of  the  DGPS  problem,  in 
particular  when  there  is  no  designated  reference  station  and  all  users  are  mobile  in  which 
case  good  relative  positioning  is  achieved.  In  addition,  treating  the  common  errors  as 
parameters  to  be  estimated  also  yields  estimates  of  the  clocks’  relative  errors.  In  relative 
GPS  good  estimates  of  the  relative  positions  of  the  rovers  are  obtained,  however,  it  is  not 
the  relative  positions,  but  rather  it’s  the  clocks’  relative  errors  that  are  directly  estimated. 

In  summary,  using  stochastic  modeling,  a  general  Kinematic  DGPS  algorithm  is 
developed  with  the  following  qualities: 

•  No  stationary  reference  station  is  envisaged  and  a  mobile  users’  formation  is 
considered.  Improved  relative  positioning  is  achieved.  This  is  often  all  one  is 
interested  in,  as  in  fonnation  flight  control. 

•  The  standard  and  network  DGPS  scenario  are  special  cases.  Improved 
absolute  positioning  is  then  achieved. 


62 


•  DGPS  is  applied  to  KGPS.  Thus,  in  this  research  DGPS  and  KGPS  are  jointly 
considered  in  a  unified  mathematical  framework. 

•  The  simulation  experiments  are  carefully  calibrated.  Obviously,  the 
positioning  accuracy  is  of  interest,  but  in  addition  the  predicted  covariance  of  the 
state/position  estimation  error  is  also  derived  and  is  compared  with  the 
experimentally  measured  estimation  error  covariance.  Having  a  correct  predicted 
estimation  error  covariance  is  particularly  important  if  GPS  measurements  are  to 
be  blended  with  measurements  provided  by  other  sensors,  as  in  INS  aiding.  The 
mathematical  result  concerning  the  probability  that  a  measurement  will  fall  inside 
the  multidimensional  l-o  error  ellipsoid  is  not  well  known,  even  in  statistical 
circles.  In  two  dimensions  it  is  always  39%,  irrespective  of  the  covariance.  This 
result  is  provided  in  the  Appendix  and  was  used  to  calibrate  the  experiments  for 
this  research. 

•  Interesting  applications  of  DGPS,  such  as  when  all  the  users  are  mobile  and 
there  is  no  stationary  reference  station  as  in  relative  DGPS,  are  listed  in  the 
conclusions. 

3. 1. 6  DGPS  Novel  Algorithms 

Algorithms  for  the  following  navigational  scenarios  are  derived: 

•  Augmented  conventional/network  DGPS  (ADGPS):  n  satellites,  one  user,  in 
(1  <  m)  stationary  reference  stations  at  surveyed  locations,  one  ( N  =  0) 
measurements  epoch. 
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•  Augmented  Kinematic  DGPS  (AKDGPS):  n  satellites,  one  user,  in  (1  <  m) 
stationary  reference  stations  at  surveyed  locations,  N  +  1  (1  <  N)  measurements 
epochs. 

•  General  Kinematic  DGPS  (GKDGPS):  n  satellites,  one  user,  m  (1  <  in)  mobile 
reference  stations,  N+  1  (1  <  AO  measurements  epochs. 

The  GKDGPS  scenario  can  also  be  viewed  as  entailing:  n  satellites,  a  team  of  in 
(, m  >  2)  mobile  users,  N  +  1  (N  >  1)  measurements  epochs.  Two  approaches  to  KDGPS 
are  developed: 

•  The  indirect  approach  to  KDGPS,  where  standard  DGPS  corrections  are 
applied  to  the  pseudoranges  and  snapshots  of  the  vehicle's  position  are  initially 
calculated.  Subsequently  the  KGPS  fitting  methodology  is  applied  to  mitigate  the 
deleterious  effects  of  measurement  noise;  and, 

•  Direct  KDGPS,  where  the  reference  receivers'  provided  pseudoranges  data 
and  the  users'  pseudorange  measurements  are  jointly  processed  in  a  centralized 
algorithm. 

Obviously,  both  in  DGPS  and  in  KGPS,  the  correct  treatment  of  common  errors  is 
most  advantageous.  This,  in  turn,  requires  a  centralized  computational  algorithm.  Linear 
regression  and  ILS  from  statistics  and  the  SVD  method  from  numerical  analysis  are  used 
with  the  payoff  being  improved  estimation  performance.  This  is  the  rationale  for  direct 
KDGPS.  In  this  research  the  emphasis  is  on  the  direct  approach. 
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KDGPS  algorithms  are  presented  where  the  ILS  approach  to  the  solution  of  the  GPS 
trilateration  equations  is  used.  The  alternative  direct  approach  [22]  to  the  solution  of  the 
GPS  trilateration  equations  can  also  be  used. 

3.2  Theory 

This  section  describes  the  theory  behind  the  development  of  the  stochastic  model- 
based  DGPS  estimation  algorithm  tested  in  the  Monte  Carlo  simulations. 

3.2.1  Overview 

Conventional  DGPS  implementations  typically  employ  a  GPS  receiver  at  a 
precisely  known  location  (reference,  or  base,  station)  to  observe  the  common  mode  errors 
in  the  pseudorange  measurements  from  each  satellite.  The  differential  corrections 
calculated  from  these  errors  are  then  transmitted  to  all  users  within  the  vicinity.  The  base 
station  is  usually  an  all-in-view  receiver  that  computes  and  transmits  corrections  for  all 
visible  satellites.  A  DGPS  remote  receiver  then  perfonns  acquisition  and  tracking  using 
nonnal  operations,  but  it  applies  the  received  corrections  to  each  pseudorange 
measurement  prior  to  computing  its  position  solution.  In  this  manner,  absolute  position 
accuracies  much  better  than  the  GPS  precise  positioning  service  (PPS)  level  are  achieved. 

The  Radio  Technical  Commission  for  Maritime  Services  (RTCM)  standard  for 
DGPS  corrections  [43]  makes  the  following  recommendations:  1)  A  reference  station 
should  not  attempt  to  remove  the  effects  of  ionospheric  and  tropospheric  errors  from  the 
broadcast  corrections.  2)  The  effect  of  reference  station  receiver  clock  error  should  be 
removed  from  the  broadcast  corrections.  3)  The  effect  of  satellite  clock  error  modeled  by 
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the  standard  polynomial  correction  applied  at  the  rover's  receiver  should  be  removed 
from  the  broadcast  corrections.  4)  The  reference  station  antenna  should  be  sited  or 
augmented  with  additional  hardware,  e.g.,  a  choke  ring,  and  the  corrections  are  processed 
to  minimize  the  effects  of  multipath. 

Also,  it  is  assumed  that  all  measurements  are  taken  simultaneously.  Latency 
issues  are  not  addressed. 

The  pseudorange  from  satellite  i  measured  by  user  k'h  receiver  is  conventionally 
modeled  as 


R 


iW  =  V(xSi  +(ySi  -yk)2+(zSi  -zj  +  Tk  +  Vjw  +  rSi  +  E, 


+  on,  +  ^trop. 


(3-2) 


where  Tk  is  the  kth  receiver  range-equivalent  clock  bias,  t_s,  is  the  i,h  satellite  range- 
equivalent  clock  bias,  E  represents  the  ephemeris  (or  SA  induced)  error,  Aion  represents 
the  ionospheric  error,  Atrop  represents  the  tropospheric  error,  and  Vi(k)represents  the 
receiver  ith  channel  measurement  error.  The  last  three  terms  in  equation  (3-2)  represent 
the  common  mode  errors.  They  are  referred  to  as  such  because  they  affect  all  receivers, 
including  the  reference  station(s),  in  a  limited  geographic  area,  in  the  same  manner. 

In  conventional  DGPS  navigation  the  reference  station  broadcasts  real-time 
(pseudorange)  corrections  that  enable  a  GPS  user  to  eliminate  the  effects  of  the  common 
mode  errors  from  his  positioning  solution.  The  extent  to  which  this  objective  is  achieved 
depends  on  the  ability  of  the  reference  station  to  separate  the  common  mode  and  non¬ 
common  mode  errors. 

A  reference  receiver  at  an  accurately  surveyed  location  (x,  y,  z)  can  calculate  the 
re fcrcnce-to-/',/'  satellite  range  as 
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(3-3) 


Ri  =A/k“x  y  +  (ySi - y  )2  +  (\~4 

The  range  differential  correction  for  the  ith  satellite  is  determined  by  differencing  the 
calculated  and  measured  reference-to-satellite  ranges: 

ADGPS  =  R,  -  Ri  _(rs,  +  Ei  +  Ai on,  +  ^tmp,  +  V i)  (3-4) 

The  sign  of  Adgps  is  motivated  by  the  RTCM  DGPS  standard  [43]  which  states  that  the 
correction  will  be  added  by  the  remote  user. 

Based  on  the  RTCM  specifications,  the  broadcast  corrections  should  be  corrected 
to  remove  the  reference  receiver  and  modeled  satellite  clock  errors.  Equation  (3-3)  shows 
the  actual  reference  station  calculation  in  conventional  DGPS  and  equation  (3-4)  shows 
the  remaining  error  sources  in  the  transmitted  correction.  Therefore,  the  broadcast 
correction  will  take  the  fonn  of  equation  (3-4).  Note  that  the  ADGps  signal  contains  the 
desired  common  mode  error  sources,  which  will  cancel  the  corresponding  errors  in  the 
user  position  calculation. 

It  is  possible  to  do  temporal  filtering  of  the  basic  DGPS  corrections  before  they 
are  broadcast  to  the  user.  The  purpose  of  the  preliminary  filtering  is  to  reduce  the  non¬ 
common  mode  receiver  measurement  noise  prior  to  broadcast.  The  basic  correction  for 
each  in  view  satellite  is  then  processed  by  its  own  Kalman  filter.  Reference  station 
algorithms  which  entail  the  temporal  filtering  of  the  basic  DGPS  corrections  before  they 
are  broadcast  to  the  user  are  not  addressed  in  this  research.  However,  in  this  research,  the 
KGPS  concept  is  synergistically  applied  to  DGPS  in  a  centralized  estimation  algorithm. 
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3.3  Novel  Kinematic  DGPS  Algorithm 


This  chapter  provides  a  detailed  description  of  the  novel  kinematic  DGPS 
algorithm  developed  to  estimate  user  position  and  velocity  for  a  network  of  users. 

3.3.1  Overview 

The  development  of  a  GPS  based  relative  position  sensing  system  is  undertaken, 
with  application  to  aircraft  or  spacecraft  formation  flying.  Accurate  real-time  solutions 
for  the  relative  vehicle  positions  are  required.  Thus,  this  research  addresses  the  general 
case  where  all  sites  are  in  motion  with  respect  to  the  earth  and  to  each  other,  a  situation 
we  refer  to  as  general  kinematic  differential  positioning  (GKDGPS).  There  are  also 
terrestrial  applications  for  which  precise  relative  positioning  is  desired  between  two  sites, 
both  in  motion  and  both  far  enough  from  fixed  land  based  sites  that  conventional  DGPS, 
as  described  in  the  previous  section,  is  not  feasible.  This  application  might  involve  the 
precise  determination  of  the  distance  between  two  ships  far  out  at  sea.  If  both  carry  GPS 
receivers,  the  precise  relative  track  between  the  two  could  be  determined  by  processing 
the  data  with  the  GKDGPS  algorithm.  The  same  applies  to  formations  of  land,  air,  or 
space  vehicles. 

Centralized  estimation  is  pursued:  The  GPS  measurements  from  all  vehicles  are 
collected  and  one  solves  for  the  state  of  the  entire  formation.  All  the  measurements  must 
be  combined  to  resolve  the  entire  fonnation  state.  Thus,  a  novel  web-based  navigation 
concept  is  developed.  This  concept  provides  optimal  estimates  of  the  formation's  state.  In 
a  large  formation,  high  communication  rates  might  be  required  in  order  to  transfer  the 
large  amounts  of  raw  GPS  data  as  is  required  in  DGPS. 
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3.3.2  Problem  Statement 


Two  or  more  mobile  users  are  considered  over  a  pre-specified,  presumably  short, 
time  interval.  The  number  of  users  and  the  specific  GPS  satellites  in  view  are  held 
constant  during  the  measurement  interval.  Only  pseudorange  measurements  are  used. 

A  specific  case  with  two  users  is  developed  to  demonstrate  the  applicability  of  the 
novel  algorithm.  At  the  beginning  of  the  time  interval,  prior  information  on  user  2's 
position  (and  velocity)  is  available.  Pseudorange  data  is  collected  over  several  epochs 
during  the  measurement  interval.  A  fixed  sampling  rate  is  used.  The  prior  information 
and  all  pseudorange  measurements  are  processed  simultaneously  at  the  end  of  the 
measurement  interval  and  estimates  of  the  position  and  velocity  of  both  users  1  and  2  are 
obtained.  The  effect  of  having  varying  degrees  of  uncertainty  in  user  2's  “reference 
station"  prior  position  information  is  examined.  In  addition,  the  effect  of  receiver  noise 
on  the  algorithm's  ability  to  correctly  estimate  the  users'  positions  and  velocities  in  the 
presence  of  user  and  satellite  clock  errors  -  the  latter,  resulting  from  Selective 
Availability  (SA)  if  it  were  turned  back  on  -  is  investigated. 

3.3.3  Theory 

The  developed  novel  algorithm  integrates  kinematic  modeling  and  differential 
GPS.  This  is  achieved  by  maximizing  the  use  of  the  information  in  the  pseudorange 
measurements  from  m  users  receiving  n  satellites  that  are  available  at  N  +  1,  relatively 
close,  time  instants  during  the  measurement  interval.  The  correct  treatment  of  the 
common  (satellite  clock)  errors  and  proper  stochastic  modeling  of  the  measurement 
situation  on  hand  are  responsible  for  achieving  improved  user  positioning. 
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A  sampling  interval  of  AT  seconds  is  used.  The  duration,  T,  of  the  measurement 
interval  consists  of  N  +  1  discrete  measurement  epochs.  Thus,  the  time  instants  when 
measurements  are  taken  are 

tj  =  AT.  j,  j  =  0, 1,...,N  (3-5) 

The  number  of  epochs  available  over  the  duration  of  the  measurement  interval  is  N  +  1 , 
where 

N  =  T  /  AT  (3-6) 

For  the  short  measurement  duration  being  stipulated,  the  users'  kinematics  are 
modeled  as  constant  speed  and  rectilinear  motion 

xk(0=  Xok  +Vk  t  (3-7) 

where  k  =  1,  2,...,  m  and  in  is  the  number  of  users;  xk,  xok,  Vk  e  R  .  The  kinematic 

model  could  include  higher  order  tenns  such  as  acceleration,  jerk,  or  more  importantly, 
could  represent  a  more  complex  motion,  e.g.,  a  satellite  trajectory  parameterized  by  its 
orbital  parameters.  However,  for  the  generic  scenario  examined  in  this  research,  equation 
(3-7)  will  be  the  kinematic  model  chosen  for  the  users’  motion. 

The  pseudorange  from  the  k'h  user  to  the  ith  satellite  at  time  t  is  modeled  as 

Pik\t)  =  <J[xk  (0  -  xSi  (of  +  [>’,  (0  -  ySi  (t)J2  +  k  (0  -  zSt  (t)f  +  Tk  +  rSi  (3-8) 
and  the  measurement  equation  is 

ZiW(t)=piW(t)  +Vi(k\t)  (3-9) 

where 

•  rA.  is  the  k,h user  range-equivalent  clock  bias  with  k  =  1,  2, ...,  m  . 
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•  t  is  the  ith  satellite  range-equivalent  clock  bias  with  i  =  1,  2, n  . 

•  (xs  ,  ys  ,  zs  )  is  the  ith  satellite  position  at  time  t  from  the  ephemeris  data  with 
i  =  1,  2, n . 

•  ( xk(t),yk(t),zk(t ))  is  the  kth  user  position  at  time  t  with  k  =  1,  2, m  . 

•  i  is  the  satellite  number  (the  number  of  satellites  in  view  is  n) . 

•  V;a  )(t)is  the  measurement  noise  in  the  ith  channel  of  user  kth  receiver  at  time 


t. 

The  measurement  noise  at  time  j  is  modeled  as 

v.JV^kN^ct1)  (3-10) 

and 

^(vijWvi 0  for  i*i',  j*j',k*k'  (3-11) 

Incorporating  the  kinematic  model  of  equation  (3-7)  for  the  users'  positions  for 
each  time  instant  j,  i.e.,  setting 

+  (K  xAr)7  (3-12) 

yk  (0->  To,  +  (KyAT)j  (3-13) 

zA-(0^z0t  +  (VkAT)j  (3-14) 

where  the  subscript  0  denotes  the  initial  position  at  time  t  =  0  (j  =  0),  allows  equation  (3- 
8)  to  be  written  as 

p-,r  =  )[*„,  +  (K  +  K,  + 

+  izot  +  ( KAT)J  ~  zSi  (0]2  +  Tk  (0  +  rSi 
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Note  that  in  equation  (3-15),  t  =jAT,j  =  0,  1,  N.  To  apply  a  kinematic  model  to  the 
users’  clock  errors,  set  rk  ( t )  =  rk  +  ik  AT  j,  j  =  0,1,. ..,7V  .  In  the  special  case  where  two 
users  (m  =  2)  and  five  (n  =  5)  satellites  are  considered,  the  parameter  of  interest  is 


6  =  [x0  , y0  , z0i ,  VXi  AT,V  AT,  Vz  AT, x0i , y 0, ,  z  V,  AT,  Vy  AT,  K  AT,  t,  , 

(3-16) 


^TSlTS2TS3Ts4Ts5TiAfT2Ar] 


where 

•  ( x0t  ,  y0t  ,  z0t  )  represents  the  k,h  user  initial  (at  time  t  =  0)  ECEF  coordinates, 
k  =  1,2,  ...,  m. 

•  ( V  ,  V  ,  V  )  represents  the  kth  user  velocity,  k  =  1,2,  . . . ,  in. 

•  rk  represents  the  kth  user  range-equivalent  receiver  clock  error,  k  =  1,2,  ..., 


in . 

•  Tk  represents  the  kth  user  velocity-equivalent  receiver  clock  drift,  k  =  1,2,  .., 
in. 

•  rs  represents  the  ith  satellite  range-equivalent  clock  error,  i=  1,2,  ...,«. 

21 

The  parameter  vector  shown  in  equation  (3-16)  is  6  e  R“  .  In  the  general  case, 
the  parameter 

6  e  r8w+"  (3-17) 

In  the  algorithm,  the  pseudoranges  received  are  composed  as  follows.  First,  the 
pseudoranges  received  by  the  kth  user  from  n  satellites  at  time  instant  j  are  used  to  form 
the  n  x  1  vector 


72 


(3-18) 


Now,  composing  the  received  pseudoranges  over  the  N  +  1  time  instants  yields 
the  n  (N  +  1)  x  1  vector 


Z{k)  = 


7  (*) 
'0 

7  (*) 


(3-19) 


and  finally,  composing  for  the  m  users  yields  the  mn(N  +  1)  x  1  “measurement  vector" 


Z  = 


r(D 

r{rn) 


(3-20) 


The  pseudorange  expression,  p(0) ,  is  composed  similarly.  Define 

A-k\0)  =  p0t  +  (VXA  T)j-xSi(t)]2+[y0k  +  ( VyAT)j-ysit)f 
+  K  +  (KA  T)j-z(t)]2 


(3-21) 


and 


a/»  =  fJKI(Q)  +  Z"k  +  tkAT  j  +  z"Si 


w 


(3-22) 


We  compose 
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fr = 


(3-23) 


Next,  composing  the  N  +  1  time  epochs  we  obtain 


fm  = 


(3-24) 


and  finally,  composing  for  the  m  users  yields  the  function  f :  R6 


%mn(N+\) 


(3-25) 


The  vector  p  is  similarly  composed,  thus  obtaining  the  function  p(6)  :  R8 
R""'(V+i).  The  nonlinear  GPS  equations  are 


Z  =  p{&)  +  W 


(3-26) 


where  the  mn(N  +  1)  equation  error  W  represents  the  composed  measurement  noise 


vector 


(3-27) 


n’N  Jmn(N+Y)xl 
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with  covariance  provided  in  equation  (3-11) 

R,  =E(WWt)=  a2Imn(N+l)  (3-28) 

Linearization  of  equation  (3-26)  with  respect  to  the  parameter  9  at  the  Ith 
iteration,  about  the  current  parameter  estimate  9{l) ,  yields 

Z  =  />09(,))  +j£  ( 9-9(n)+W  (3-29) 

Ou  ^(/) 

Equation  (3-29)  is  rearranged  to  yield  the  linear  regression  in  9 

Z  +  0W  -  p{9U))  =  9  +  W  (3-30) 

Ou  q(1)  Ou  q(1) 

Note  that  as  a  result  of  our  carefully  performed  linearization  of  equation  (3-20),  the  LHS 
of  equation  (3-30)  consists  not  only  of  the  measurement  vector  Z,  as  is  the  case  in 
conventional  ILS,  and  two  additional  terms  are  included.  Obviously,  when  p(9)  is  linear, 
there  is  no  difference  between  our  improved  ILS  algorithm  and  the  conventional  ILS 
algorithm,  however  when  the  nonlinearity  is  strong,  the  new  ILS  algorithm  will 
outperform  the  conventional  ILS  algorithm. 

The  calculation  of  the  regressor  matrix  requires  the  composition  of  the  partials  of 
p(9)  with  respect  to  the  parameter  vector  9 ,  viz.,  let 

=  3Ai'  ’(g)  (3-31) 

o9 

where//;  jU)  is  a  1  x  (7m  +  n)  row  vector.  Its  entries  are  explicitly  given  by 


75 


where 


(xs  ,  ys  ,  zs  )  is  the  i  satellite  ephemeris  at  time  j. 


2 

eik,el  ( j),e2  (y),e3  (j)  are  8m+  4n  row  vectors  of  zeroes,  with  l's,  j,  and j 


located  at  positions  indicated  in  their  subscripts,  according  to 


e\ik  ( J )  —  eTk-6Jk-3,lm+4i-3 
e2 i  k  O')  —  ^7 k—5,1  k—2,1  m+Ai-2 


ei,k  eik,lm+Ai 


•  The  matrix 

0,k  0)  =  %  r  (./)  e,.t  (./)  +  e2.  t  r  (j)e2lk  ( j )  +  %  r  0>3a  0)  (3-33) 

The  composition  of  the  regressor  matrix  is  similar  to  the  process  employed  for  Z  and 
f.  First,  for  the  n  satellites,  the  n  x  ( 7m  +  n)  matrix 

ix/»i 


H^k\d)= 


(3-34) 


Hjk\e ) 


is  fonned,  followed  by  composition  over  the  N  +  1  time  epochs  yielding  the  n(N  +  1)  x 
(7m  +  n)  matrix 


H*\0) 


H(k)  (6)  = 


(3-35) 


HN(k\6) 
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Next,  set 


H(k\0 )  :=  [. H(k\0 )  |  //,] 

where  all  the  entries  of  the  n(N  +  1)  x  n  matrix  Hi  are  0,  except  its  kth  column  which  is  (0, 
1,2,  . . .,  AO  T,  and  finally,  for  the  m  users,  the  mn(N  +  1)  x  (8m  +  n)  regressor  matrix 

Hm(0) 

H(0)=  •  (3-36) 

H(m)(0 ) 

is  obtained.  Thus  Equation  (3-29)  yields  the  linear  regression 

z  +  e(1)  -  f(0(I) )  =  -^  o  +  w 

SO  dO  yo  (3-37) 

=  H(60))6  +  IT 

The  function  p(O) :  r8",+  ''  — *  R""'(  V  1 1  is  linear  in  the  users'  and  satellites'  clock  error 
parameters  and  therefore  the  function  p(0)  is  replaced  by  the  function  / (6) in  the  LHS 
of  equation  (3-37). 

The  solution  for  the  users'  “positions",  viz.,  the  parameter  vector  0 ,  may  be  obtained 
using  least  squares,  whereupon  the  Iterated  Least  Squares  (ILS)  algorithm  is  established: 

0(M]  =  [HT  (0W)  H(0(l))Yl  HT  (0(ly)[Z  + 

dO 

where 


0(I)  ~  f0U))]  (3-38) 

M) 


•  Z  is  the  received  “stacked  up"  pseudorange  measurements  mn(N  +  1)  x  1 
vector,  according  to  equation  (3-20). 

•  f  is  the  mn(N  +  1)  x  1  vector  of  ranges  calculated  according  to  equation  (3-21 ). 


77 


•  H  is  the  mn(N  +  1)  x  (7m  +  n)  regressor  matrix  formed  from  the  partials  of  p 
according  to  equation  (3-36). 

3.3.4  Reduced  parameter  vector 

The  regressor  H  is  further  examined.  As  stated  in  equation  (3-16),  the  parameter 
vector  contains  (8m  +  n)  variables:  3  position  components,  3  velocity  components,  a 
clock  bias  variable,  and  a  clock  drift  variable  for  each  of  the  two  (m  =  2)  users,  as  well  as 
the  five  (n  =  5)  satellite  clock  errors  -  21  variables  in  total.  The  main  objective  is  to 
estimate  the  (non-reference)  users'  position  and  velocity  rather  than  the  users'  and 
satellites'  clock  errors.  With  this  in  mind,  the  regressor's  matrix  structure  is  examined  and 
the  algorithm  is  modified  according  to  the  following  analysis.  To  simplify  the  exposition, 
the  users  clock  drifts  are  not  considered  herein  and  the  reduced  parameter  vector  has  the 
dimension  7m  +  n 

Define  the  mn(N  +  1)  x  (m  +  n )  matrix 

B  =  [bl,b2,bi,...,cl,c2,ci,...,cn]  (3-39) 

where 

•  bk  is  the  column  of  H  operating  on  user  klh  clock  error  variable,  k  =  1 ,  2, . . . ,  m. 

•  ck  is  the  column  of  H  operating  on  satellite  i,h  clock  error  variable,  i=  1,2,..., 
n. 

Indeed,  the  structure  of  the  regressor  matrix  H(6{,) )  is 

H  =  [Diag(Hm0V)yt U  I B]  (3-40) 
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where  the  n(N  +  1)  x  6  matrices  H(k)  are  the  H(k)  matrices  with  only  the  six  columns 
(6k  -  5),  (6k  -  4),  (6k  -  3),  (6k  -  2),  (6k  -  1)  and  6k,  retained;  e.g.,  in  the  case  where  m  = 
2,  and  therefore  k  =  1,2 


H  = 


' Ha) 

vO 


0 

h{2) 


B 


(3-41) 


where  the  H 1 1 1  and  // 1 2 1  matrices  are  explicitly  given  by  the  partials  of  f  in 
(Xo^Zo^A^AT^AT)  and  (x0i,y02,z02,VXiAT,VyiAT,VZ2AT) ,  respectively, 


and  the  matrix  B  is  given  by 


B  = 


D^g(en(N+l))mk=t 


where 

•  en(N+l)  is  a  vector  of  ones  of  size  (N  +  1)  x  1  . 

•  I n  is  an  identity  matrix  of  size  n. 

For  example,  for  m  =  2, 


B  = 


cn(V+ 1) 

0 


0  ' 

en(N+ l)y 


(3-42) 


(3-43) 


The  following  holds 
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(3-44) 


Yjbk  =  Yjci 

k= 1  ;=1 

Obviously,  the  matrix  B  (and  therefore,  the  regressor  H)  is  rank  deficient.  The  rank 
deficiency  is  exactly  1.  Thus,  perfonn  the  full  rank  factorization  of  B, 

B  =BXK  (3-45) 

where  Bx  is  full  rank  (m  +  n  -  1)  and  is  an  mn(N  +  1)  x  (m  +  n  -  1)  matrix  of  the  form 


B\  \bx ,  b2 .,  bm_x ,  ^  bk,cx,c2,...,  cnX  ] 

k= 1 


(3-46) 


For  example,  in  the  special  case  where  m  =  2 


B{  =[bl,bl+b2,cl,c2,...,cn_l\  (3-47) 

This  choice  of  basis  is  responsible  for  inserting  a  column  of  ones  in  the  regressor  matrix 
H.  In  the  parlance  of  linear  regression,  an  “intercept"  variable  is  then  included.  The  latter 
has  the  beneficial  effect  of  absorbing  truncation  errors  caused  by  the  linearization  of  the 
RHS  of  equation  (3-30).  This  basis  choice  also  has  the  effect  of  yielding  the  estimates  of 
clock  error  differences,  as  indicated  in  equation  (3-37);  however,  we  are  mainly 
concerned  with  the  estimation  of  user  l's  position  and  velocity  parameters. 

Solving  equation  (3-45)  for  K  yields  the  blocked  (m  +  n  -  1)  x  (m  +  n)  matrix 


K  = 


I 

0 


n 


( n—\)xm 


0 


mx(n- 1) 


/ 


n- 1 


(3-48) 


where 


/ 


n 


I 

0 


m- 1 


b  = 


0 


(m-l)xl 


1 


(3-49) 

(3-50) 
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and 


•  0  is  a  zeroes  matrix. 

•  en  is  a  vector  of  ones  of  length  n. 

Next,  partition  the  regressor 

H  =  [H\B]  (3-51) 

where  H  =  Diag(IIU]  (0'i]  ))"'=l  is  an  mn(N  +  1)  x  6m  matrix  consisting  of  the  columns  of 

H  operating  on  the  users'  position  and  velocity  parameters  only.  Also  define  the  reduced, 
full  rank,  matrix 

(3-52) 

Next,  perform  the  full  rank  factorization  of  the  regressor  H  (=  \H  B]) 

H  =  H{KX  (3-53) 

i.e.,  the  following  equation  is  solved  for  the  (7m  +  n  -  1)  x  (7m  +  n )  matrix  K]  : 

[h\b,K]=  [h\b,]K,  (3-54) 

We  calculate: 


±  6m  w  6mx(m+n) 

0  K 

(m+n—\)x6m  ( m+n-l)x(m+n ) 


(3-55) 


Finally,  the  reduced  parameter  vector  is  defined 

6>  =KX0  (3-56) 

The  reduced  parameter  vector  0X  e  7"'+""'  consists  of  the  users'  position  and  velocity 
parameters  we  are  interested  in,  as  well  as  linear  combinations  of  the  user  and  satellite 
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clock  errors.  For  the  specific  scenario  examined  in  which  we  have  two  mobile  users  and 
five  satellites  visible,  this  yields  the  (18  x  1)  parameter  vector 


Ox  =  K ,  To, ,  ,  K,  A T,  vyi  AT ,  Vz<  AT ,  x^  ,y^,z0,VxAT,VyAT,VzAT 


Tx  —  X  o  «,  X  ^  +  X  *  X  —  X  +  X  —  X  ,  X  —  X  ,  X  —  X  1 

1  2^2  s5  ’  ^  s2  's5’  's3  85^84  s<  J 


(3-57) 


Since 


H6  =  HXKX6  =  HX9X  (3-58) 

equation  (3-37)  is  written  with  the  reduced  parameter  6X  : 

Z  +  H2(02O))02il)  -  f{9{l))  =  Hx(92l))9xl)  +  W  (3-59) 

where,  in  addition,  the  further  reduced,  (12  x  1),  parameter  vector  is  used 


ei  =  [-S » Zo, » zo,  ’K,A  T,  Vyi  AT,VzAT,x02,y0i,z02,VXiAT,  Vyi  AT, 
VZA  T]t 


(3-60) 


Thus,  the  further  reduced  parameter  vector  02  is  stripped  of  the  clock  error 
parameters  of  0X  ,  and  02  (and  not  9X )  is  used  on  the  LHS  of  equation  (3-59)  because  the 
function  f  is  not  dependent  on  the  time  parameters.  Accordingly,  the  matrix  H2  is 
composed  of  the  first  6m  columns  of  H  associated  with  the  parameters  featuring  in  A, 
(positions  and  velocities,  no  clock  errors).  Thus,  H2  =  Diag(lTk>)"=l  . 

With  reference  to  equation  (3-57)  it  is  noted  that  not  using  differencing  and 
double  differencing  as  is  standard  practice  in  DGPS  and  instead  treating  the  common 
errors  as  parameters  to  be  estimated  yields  estimates  of  the  clocks’  relative  errors.  In 
relative  DGPS  good  relative  position  estimates  are  obtained  but  interestingly,  it  is  not  the 
relative  positions,  but  rather,  the  clocks’  relative  errors  that  are  directly  estimated.  Should 
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laser  ranging  be  available,  as  in  the  European  GIOVE-A  experiment,  clock  and 
ephemeris  errors  could  be  separated. 

3.3.5  Prior  Information 

The  linear  regression  featuring  in  the  ILS  algorithm  can  be  augmented  to  include 
prior  information  on  user  2  (position  and  velocity),  as  in  the  conventional  and  network 
DGPS  [36,  38,  42].  The  prior  information  is  provided  in  the  form 

*o2  e  N(x0i  >^2*o2) 

Vo2  e  N(y02  ,a\02) 
zo,  e  N(zo2 

-  ,  (3-61) 

K2  e  N{VXi  ,<J2v„ ) 

VV2  e  N(Vyi  ,<j\2) 

VZl  e  N(VZi  ,a\2) 


Using  very  large  a  parameters  in  equation  (3-61)  is  tantamount  to  the  stipulation  that  no 
prior  information  on  user  2  initial  state  is  available. 

The  linear  regression  in  equation  (3-59)  is  now  augmented  as  follows: 


Z  := 


\Z\) 


(3-62) 


where 


Zi  = 


fx  ^ 
A02 

To, 

*o2 

K 

Vy* 

^  J 


(3-63) 
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In  addition,  the  regressor  H1  is  augmented 


H,  :=  — 
M 


(3-64) 


where  the  6  x  (7m  +  n  -  1)  selector  matrix  is,  e.g.,  in  the  case  where  m  =  2, 


M  =  [0  |0 


6x(m+n- 1)  - 


(3-65) 


and  when  appended  to  //,  picks  out  the  x,t  ,  v0  i ,  z()  ,VxAT,VyAT,V,AT  elements  of  the 


parameter  vector  0X  .  Moreover, 


H2  :=  JA 
0 

n' 


(3-66) 


(3-67) 


m«(A^+l)+l 


mn(N+\)+6 


(3-68) 


Additionally,  a  weighting  matrix  R  is  included  to  correctly  incorporate  the 
confidence  level  in  the  “reference  "  receiver's  (user  2)  prior  information  on  position, 
velocity,  and  possibly,  range  equivalent  clock  error: 


R  =  Diag(Rl,R2 ) 


(3-69) 


where 


-^1  I  mn(N+\) 


(3-70) 
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is  determined  by  the  measurement  noise  variance  cr,  and  the  diagonal  matrix  R2 

contains  the  prior  infonnation  data,  viz.,  the  standard  deviations  of  the  reference  station's 
initial  position  (and  velocity) 

a2x02  0  0  0  0  0 

0  cr2j02  0  0  0  0 

0  0  cr2z0,  0  0  0 

0  0  0  rx2  0  0  (3-71) 

0  0  0  0  a2vy2  0 

0  0  0  0  0  cx2vzi 


where 

o-2X02  =  E(W 2  mn(N+l)+\ ) 

CT  =  E(JV  mn(N+\)+2 ) 
o-2-v2  =  E(W 2  mn(N+ 1)+3  ) 
cr\2  =  E(W 2  mn(N+\)+4  ) 
cr2r„  =  E(W2  mn(N+\)+5  ) 

<r\  =  E(W2  mn(N+\)+6  ) 

This  finally  yields  the  enhanced  ILS  algorithm 

=  [//  7 {H'n)R  7/l(«  ;,,|]  'll  riH"!lR  [Z 

+  H20m)6m  -  /(O'1')] 

and 


1  =  0,  1, 


(3-72) 


(3-73) 


Finally,  if  the  satellite  clock  errors  are  known  to  high  precision,  as  provided  by 
NASA’s  Internet-Based  Global  Differential  GPS,  and  the  residual  errors  can  be  isolated 
to  the  ionospheric  delay,  then  the  additional  prior  infonnation  Ajon  «7V( 25, 225)  can  be 
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used  based  on  the  ionospheric  retardation.  That  is,  the  linear  regression  is  augmented  to 
include  the  equations 

25  =  Aionj  +vp  vf  *  N(0,  255),  i  =  l,...,n 

3.4  Simulation  Results 

This  section  presents  the  simulation  results  and  validates  the  algorithm's  ability  to 
correctly  estimate  the  parameters  of  interest. 

3.  4.1  Simulation  Scenarios 

In  all  the  simulation  scenarios  two  (m  =  2)  mobile  users  flying  along  parallel 
tracks,  approximately  180  nautical  miles  (LEO)  above  the  earth  surface  and  separated  by 
10,000  meters,  are  considered.  The  number  of  satellites  in  view  is  n  =  5  and  n  =  8  to 
illustrate  the  effect  of  satellite  availability.  The  normalized  user  velocities  V  AT  are 
VY  =  Vr  =  100m /sec  and V„  =  Vt  =V_  =  V.  =0.  In  the  simulation  experiments 

the  users'  range-equivalent  clock  bias  are  fixed  at  t1  =  r,  =  200m ,  and  all  the  satellites 
are  subject  to  the  same  range-equivalent  clock  errors:  rs  =  100m,  /  =  1,2,...,  n  or 
rs  =  30m,  i  =  1,2,...,// .  The  number  of  measurement  epochs  considered  varied  between  N 

=  1  and  N  =  11.  In  the  case  where  N  =  1 ,  kinematic  GPS  is  not  possible  and  the  users' 
velocities  cannot  be  estimated.  This  scenario  does  however  correspond  to  the 
conventional  stand-alone  DGPS  paradigm  and  allows  for  a  direct  comparison  to  be  made 
of  the  perfonnance  of  the  novel  algorithm  and  the  differencing-based  conventional  DGPS 
algorithm. 
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The  random  residual  errors  and  receiver  noise  levels  used  in  the  simulation  have  a 


of  0.75m,  1.50m,  or  5.00m.  Ten  iterations  (T=10)  for  the  ILS  algorithm  were  performed 
for  all  scenarios. 

In  the  simulation  experiments,  the  actual  user  2  initial  position  is  randomly 
chosen  according  to  the  specified  normal  distribution;  the  latter  gauges  our  confidence  in 
the  user  2  prior  position  information.  Both  the  nominal  user  2’s  initial  position  and  the  a 
information  are  provided  to  the  DGPS  algorithm,  but  not  the  realized  (actual)  user  2 
position  used  to  generate  the  pseudorange  data  given  to  the  algorithm.  The  algorithm's 
estimation  performance  is  gauged  by  measuring  the  distance  of  the  computed  position 
estimate  from  the  actual  true  position  which,  of  course,  is  known  during  the  controlled 
simulation  experiment.  This  applies  to  both  user  1  and  user  2.  For  each  simulation 
scenario,  100  Monte  Carlo  (MC)  runs  were  performed. 

3.4.2  Numerical  Results 

The  DGPS  algorithm's  ability  to  accurately  estimate  the  position  of  user  1  in  the 
face  of  large  satellite  clock  errors  is  strongly  dependent  on  the  accuracy  of  the  reference 
station's  position  information.  This  is  in  particular  true  when  the  number  of  satellites  in 
view,  n,  is  5.  When  the  number  of  satellites  in  view  is  increased  to  n  =  8,  the  estimation 
accuracy  improves  considerably,  even  in  the  case  of  one  measurement  epoch  (N  =  1) 
only.  This  is  to  be  expected:  when  N=  1,  the  number  of  unknowns  is  4m  +  n  1=4x2  + 
n  1  =  7  +  n.  The  number  of  equations  is  mn  =  2n.  Evidently,  to  obtain  a  navigation 
solution  using  the  novel  algorithm  we  need 

7  +  n  <  2n  (3-74) 

i.e., 
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7  <n 


(3-75) 


At  the  completion  of  the  MC  experiment,  the  kt  l  user  position  estimation  error  is: 


100  y 


100 

2 

100 

2 

100 

ZW  -*of) 
p=l 

+ 

5X  -Tof) 

p=i 

+ 

IX -4f) 

p= i 

(3-76) 


where,  unbeknown  to  the  estimation  algorithm,  (jc0  ,  y0k ,  z0  )  is  the  true  initial  position  of 
user  k  and  is  used  to  generate  the  data  in  the  simulation,  which  the  estimation  algorithm 
operates  on,  and  (x^\ 's  the  algorithm-produced  position  estimate  in  the  pth 

MC  run.  The  estimation  error  spread  is  gauged  by  the  estimation  error  variance.  The 
latter  is  provided  by  the  estimation  algorithm  and  is  referred  to  as  the  predicted 
estimation  error  variance  cr  ,k  =  1,2  : 

Pk 


aPi  =  yj((Hx  R  XHX  r%  +  ((HTxRlHx  +  (3-77) 

and  similarly 

=^((H{R-XHX  r‘)7  7  +((HfRlHi  y\&+((H^RxHx  Y\g  (3-78) 


Finally,  the  experimentally  obtained  estimation  error  variance  is  calculated  as 


1/100  ^  100  ^  1UU 

+Z(Tot->'of))  +Z(zot-fof))  (3-19) 


p=i 


p= i 


p= i 


The  performance  of  the  novel  estimation  algorithm  is  summarized  in  Tables  4  and  5  on 
the  next  page.  The  user  1  position  estimation  error  is  less  than  4  m  even  though  the  two 
users'  range  equivalent  clock  errors  are  rx  =  r2  =  5m ,  the  eight  satellites'  large  range 
equivalent  clock  errors  are  rs  =  1 00 in,  i  =  1,2,..., 8,  the  inaccuracy  in  the  reference  station 
(user  2)  position  is  5m  and  the  relatively  large  receiver  noise  variance  is  5m. 


88 


Table  4.  Novel  Algorithm  Estimation  Error  &  Accuracy  [meters]  with  m  =  2 ,n  = 

8,  N=  1,  L  =  10  for  1  epoch 


Referenced  o 

Noise  o 

ei 

Opi 

Oei 

0.0001 

0.75 

0.501 

1.583 

2.874 

0.0001 

1.50 

0.875 

2.945 

3.459 

0.0001 

5.00 

0.841 

9.599 

8.890 

0.75 

0.75 

0.502 

2.003 

5.204 

0.75 

1.50 

1.571 

3.165 

5.642 

0.75 

5.00 

0.353 

9.664 

10.489 

1.50 

0.75 

0.907 

5.151 

16.160 

1.50 

1.50 

2.944 

5.615 

16.442 

1.50 

5.00 

1.558 

10.550 

20.001 

5.00 

0.75 

2.535 

10.073 

31.425 

5.00 

1.50 

2.913 

10.302 

33.069 

5.00 

5.00 

3.930 

13.354 

36.874 

The  estimation  error  is  significantly  reduced  when  five  measurement  epochs  are 
used.  In  this  case,  and  when  the  uncertainty  in  the  user  2  position  is  small  (Ref.  cr  = 
0:0001),  the  following  result  is  obtained  in  Table  5: 
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Table  5.  Novel  Algorithm  Estimation  Error  &  Accuracy  [meters]  with  m  =  2 ,n  = 

8,  N=  5,  L  =  10  for  5  epochs 


Referenced  o 

Noise  o 

ei 

Opi 

Oei 

0.0001 

0.75 

0.225 

0.523 

0.714 

0.0001 

1.50 

0.255 

1.743 

1.508 

0.0001 

5.00 

0.285 

3.486 

5.129 

Finally,  the  conventional  DGPS  algorithm's  results  are  given  in  Table  6.  As 
expected,  the  estimates  in  Table  5  are  better  than  the  estimates  in  Table  6.  This  applies  to 
both  the  estimation  error,  and  to  the  estimation  error  variance. 


Table  6.  Conventional  DGPS  Estimation  Error  &  Accuracy  [meters]  with  m  =  2,  n 

=  8,  N  =  1,  L  =  10  for  5  epochs 


Referenced  o 

Noise  o 

ei 

Opi 

Oei 

0.0001 

0.75 

0.290 

1.355 

0.983 

0.0001 

1.50 

0.355 

1.560 

2.172 

0.0001 

5.00 

0.528 

2.486 

7.100 

Furthermore,  with  n  =  8  satellites  in  view  it  is  possible  to  obtain  a  good  estimate  of  the 
users'  relative  positions  even  without  prior  information  on  the  reference  station's  (user  2) 
position  -  as  required  in  formation  control. 
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IV  Conclusions 


This  research  presents  the  theories,  models,  and  simulation  results  for  a  direct, 
closed-fonn  and  efficient  new  position  determination  algorithm  developed  for  a  stand¬ 
alone  user  and  the  stochastic  model-based  DGPS  estimation  algorithm  constructed  for  a 
network  of  mobile  users.  The  theoretical  developments  stochastically  modeled  the 
pseudorange  measurement  with  random  white  noise  and  solved  for  position  as  a 
stochastic  estimation  problem.  The  following  sections  present  the  conclusions  obtained 
from  this  research. 


4.1  A  Direct,  Closed-Form  and  Efficient  New  Position 
Determination  Algorithm 

The  performance  of  the  novel  two-step  algorithm  is  comparable  to  that  of  the 
baseline  ILS  algorithm.  Furthermore,  it  retains  all  the  attractive  features  that  motivated 
the  development  of  the  closed-form  algorithm  in  the  first  place.  Hence,  considering  the 
closed-fonn  algorithm  supplemented  by  the  Kalman  update  algorithm  as  a  single  two- 
step  GPS  position  determination  algorithm,  a  novel  algorithm  with  the  following 
attributes  has  been  developed: 

•  The  performance  under  conventional  navigation  scenarios  is  comparable  to 
that  achieved  by  the  conventional  ILS  algorithm. 
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•  The  algorithm  is  self-contained;  hence  it  can  be  used  under  any  geometrical 
conditions  without  the  need  for  externally  provided  initialization,  and  a  degree  of 
autonomy  is  thus  achieved. 

•  The  algorithm  is  computationally  efficient  because  it  requires  less  iteration 
than  the  ILS,  and  its  FLOPS  count  is  lower. 

•  The  algorithm  has  the  capability  to  produce  a  data-driven  estimate  of  the 
receiver  measurement  noise  strength  ( a  )  when  the  number  of  satellites  in  view  n 
>  6  and  is  able  to  predict  its  estimation  error  covariance. 

•  Good  estimated  performance  under  conditions  of  high  GDOP  is  obtained. 
The  horizontal  positioning  performance  of  the  novel  two-step  algorithm  under 
poor  geometry  conditions,  e.g.  when  ground-based  planar  arrays  of  pseudolites 
are  used,  is  similar  to  that  of  the  conventional  ILS  algorithm.  Moreover,  there  are 
no  restrictions  on  the  user  position  and  an  initial  user  position  guess  is  not 
required  as  long  as  the  user  was  within  the  confines  of  the  outer  radius  of  the 
circular  pattern.  This  may  prove  beneficial  in  testing  range  applications  where  the 
conventional  iterative  algorithm  is  at  risk  of  failure,  imposing  restrictions  on  the 
flight  test  trajectory  and  altitude. 

In  conclusion,  the  benefits  of  the  two-step  algorithms  are  computational 
efficiency,  data-driven  predictions  of  the  noise  strength  of  the  pseudorange  measurements 
and  the  estimation  error  covariance,  no  need  for  an  initial  position  guess,  and  may 
provide  better  performance  under  poor  geometry. 
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4.2  Stochastic  Model-Based  DGPS  Estimation  Algorithm 


In  this  dissertation  a  network-based  navigation  concept  is  advanced  and  a  novel 
algorithm  for  KDGPS  data  processing  is  presented.  Specifically,  the  accurate  relative 
(and  absolute)  positioning  of  a  team  (or  fonnation)  of  mobile  vehicles  is  considered.  The 
measurement  scenario  is  correctly  modeled,  that  is,  a  stochastic  framework  is  developed 
where  rather  than  using  double  differencing  the  common  errors  are  explicitly 
acknowledged  and  a  novel  centralized  estimation  algorithm  is  rigorously  derived.  The 
novel  DGPS  algorithm  estimates  the  users'  positions  and  velocities  utilizing  only 
pseudorange  measurements.  One  of  the  benefits  of  the  novel  algorithm  is  its  ability  to 
account  for  the  uncertainty  in  the  reference  receiver's  position  and  allows  for  the 
inclusion  of  the  information  on  the  reference  user's  position  accuracy.  Indeed,  DKGPS 
exploits  the  same  self  calibration  concept  as  DGPS  to  provide  comparable  accuracies  in  a 
dynamic  environment.  DGPS  can  provide  enhanced  positioning  accuracies  (on  the  order 
of  lm)  by  using  a  stationary  base  station,  a.k.a,  reference  station,  at  a  surveyed  location. 
Unlike  DGPS,  KDGPS  does  not  require  a  designated  base  station  and  consequently 
cannot  improve  absolute  accuracies.  Instead,  it  provides  a  similar  degree  of  relative 
position  accuracy  between  the  mobile  GPS-equipped  platforms  as  required  in  fonnation 
control.  A  data  driven  estimate  of  the  predicted  covariance  of  the  position  estimation 
error  is  provided.  Last  but  not  least,  insight  into  the  structure  of  the  estimation  problem 
on  hand  is  gained.  In  relative  DGPS  good  estimates  of  the  relative  positions  of  the  rovers 
are  obtained;  however,  it  is  not  the  relative  positions,  but  rather,  it’s  the  clocks’  relative 
errors  that  are  directly  estimated. 
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Simulations  entailing  MC  experiment  validate  the  algorithm's  ability  to  estimate 
the  variables  of  interest.  The  results  are  also  compared  to  conventional  DGPS  scenarios, 
where  the  reference  station(s)  is  (are)  stationary  and/or  only  one  set  of  pseudorange 
measurements  is  available.  The  positioning  accuracy  improvements  achieved  with  the 
novel  algorithm  are  gauged  against  the  performance  of  conventional  DGPS. 

The  generic  nature  of  the  mathematical  development  presented  in  this  dissertation 
which  entails  stochastic  modeling,  the  use  of  the  method  of  linear  regression  from 
statistics,  and  the  adaptation  of  the  SVD  method  from  numerical  analysis  are  conducive 
to  a  broad  range  of  applications  of  the  estimation  algorithm.  These  entail 

1)  Augmented  conventional  DGPS. 

a.  Correctly  treat  measurement  noise  effects  in  conventional  DGPS,  where 
differencing  causes  the  equation  error  variance  to  increase  to  2  cr 2  and  correlation  is 
introduced.  Furthermore,  the  equation  error  is  no  longer  homoscedastic  and  a  weighted 
ILS  must  be  used. 

b.  Treat  ephemeris  errors  (not  just  satellite  clock  errors)  as  unknown  parameters 
to  be  estimated. 

2)  Formations  of  multiple  spacecraft  are  envisioned  for  many  planned  space 
missions,  and  GPS  will  play  an  important  role  as  a  sensor.  Tight  control  of  the  vehicles’ 
positions  relative  to  each  other  is  very  important.  Formation  technologies  for  spacecraft 
will  enable  multiple  distributed  spacecraft  to  act  in  a  unified  manner.  This  will  enable 
new  scientific  missions  involving  distributed  but  coordinated  measurements,  leading  to 
improved  stellar  interferometry,  gravimetry,  and  synthetic  aperture  radars.  GPS,  and,  in 
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particular  generalized  KDGPS,  provides  a  promising  relative  navigation  sensor.  Similar 
success  is  foreseen  in  aircraft  fonnation  flying  applications. 

3)  Relative  GPS  guidance,  where  a  targeting  system  provides  a  moving  target’s 
coordinates  using  an  on-board  GPS  receiver  for  aiding  its  inertial  reference  system.  A 
GPS  receiver  is  also  on  the  weapon  to  track  the  same  satellites  as  the  GPS  receiver  on  the 
targeting  aircraft.  Thus  some  of  the  GPS  bias  errors  in  the  targeting  data  are  removed  and 
guidance  accuracies  can  be  improved  to  near  “precision"  level  without  the  use  of  a 
tenninal  guidance  seeker. 
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Appendix  A:  Acronym  List 


AFIT 

Air  Force  Institute  of  Technology 

AFRL 

Air  Force  Research  Laboratory 

CPGPS 

Carrier  Phase  Global  Positioning  System 

DCM 

Direction  Cosine  Matrix 

ECEF 

Earth-Centered,  Earth-Fixed 

EKF 

Extended  Kalman  Filter 

GPS 

Global  Positioning  System 

INS 

Inertial  Navigation  System 

PR 

Pseudorange 

RMS 

Root-Mean-  S  quare 

SA 

Selective  Availability 

sv 

Satellite  Vehicle 

GDOP 

Geometric  Dilution  Of  Precision 

UERE 

User  Equivalent  Range  Error 
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Appendix  B:  Bivariate  Gaussian  Distribution 


To  evaluate  the  performance  of  the  novel  algorithm,  the  probability  of  the  position 
estimate  to  be  within  an  ellipsoid  was  computed.  For  a  Bivariate  Gaussian  Distribution, 
the  probability  of  the  position  estimates  within  the  la  can  be  obtained  as  follows: 
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2^-|det(f5) 


11 


~-xTPx 
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dx 


2  X=  x  \x‘  P~'  x  <  1 


First,  let  the  estimate  error  covariance  matrix  to  be 


P  =  T 1 


A,2  0 
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where  T  is  the  orthonormal  transformation  matrix  for  P.  The  inverse  of  the  error 
covariance  matrix  can  be  written  as: 
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then  obviously 
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Thus,  the  transformation  from  y  into  x  is 
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where 


det(r)  =  At  X, 

Using  the  above  relationships,  the  probability  of  the  position  estimate  to  be  within  the 
ellipsoid  can  be  found  as 
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